From ae2b6fe4ae898aa90705e3f13ea3aeec716f77bb Mon Sep 17 00:00:00 2001 From: Nanda Date: Wed, 4 Nov 2020 12:40:35 +0700 Subject: [PATCH 1/5] :construction: correct some suggestions --- include/cell.tcc | 2 +- include/node.tcc | 21 ++++++++++----------- 2 files changed, 11 insertions(+), 12 deletions(-) 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/node.tcc b/include/node.tcc index e4840f1a0..d819cd9f2 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()); From b8a7d89ccd7c6c7756e57494c6467850f4b67a7d Mon Sep 17 00:00:00 2001 From: Nanda Date: Wed, 4 Nov 2020 16:28:45 +0700 Subject: [PATCH 2/5] :construction: move multiphase functions outside mesh.tcc --- include/mesh.h | 27 +- include/mesh.tcc | 495 ------------------------------------ include/mesh_multiphase.tcc | 494 +++++++++++++++++++++++++++++++++++ 3 files changed, 515 insertions(+), 501 deletions(-) create mode 100644 include/mesh_multiphase.tcc 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 From ae92ec076608a369c3752aae10887aacb2747de7 Mon Sep 17 00:00:00 2001 From: Nanda Date: Wed, 4 Nov 2020 16:29:03 +0700 Subject: [PATCH 3/5] :construction: add groups in cell and node functions --- include/cell.h | 28 ++++++++++++++++++++++------ include/node.h | 13 +++++++++++-- 2 files changed, 33 insertions(+), 8 deletions(-) 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/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 From 62e065ccab8bdf786b87db6bb78e867774857baf Mon Sep 17 00:00:00 2001 From: Nanda Date: Wed, 4 Nov 2020 16:29:22 +0700 Subject: [PATCH 4/5] :dart: modify test to fix error in assertion --- tests/mesh_test_2d.cc | 7 ------- tests/mesh_test_3d.cc | 7 ------- 2 files changed, 14 deletions(-) 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 From 32f1c5a98d5689f54730d70a6b3e5d9006637321 Mon Sep 17 00:00:00 2001 From: Nanda Date: Thu, 5 Nov 2020 00:03:06 +0700 Subject: [PATCH 5/5] :construction: :dart: change nPhase to NPhase in nodes --- include/node.tcc | 2 +- include/node_base.h | 8 +- include/node_twophase.tcc | 92 +- include/solvers/mpm_explicit_twophase.tcc | 32 +- tests/nodes/node_test.cc | 18 +- tests/nodes/node_twophase_test.cc | 996 +++++++++++----------- tests/particles/particle_twophase_test.cc | 70 +- 7 files changed, 609 insertions(+), 609 deletions(-) diff --git a/include/node.tcc b/include/node.tcc index d819cd9f2..2fd2bb162 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -625,7 +625,7 @@ void mpm::Node momentum = property_handle_->property("momenta", prop_id_, *mitr, Tdim); const Eigen::Matrix change_in_momenta = - velocity_.col(mpm::NodePhase::nSolid) * mass - momentum; + velocity_.col(mpm::NodePhase::NSolid) * mass - momentum; property_handle_->update_property("change_in_momenta", prop_id_, *mitr, change_in_momenta, Tdim); } diff --git a/include/node_base.h b/include/node_base.h index baae25e1a..b8e206835 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -19,10 +19,10 @@ namespace mpm { //! Particle phases enum NodePhase : unsigned int { - nSolid = 0, - nLiquid = 1, - nGas = 2, - nMixture = 0 + NSolid = 0, + NLiquid = 1, + NGas = 2, + NMixture = 0 }; //! NodeBase base class for nodes diff --git a/include/node_twophase.tcc b/include/node_twophase.tcc index 3d153a0eb..d1f6f033e 100644 --- a/include/node_twophase.tcc +++ b/include/node_twophase.tcc @@ -19,26 +19,26 @@ bool mpm::Node:: compute_acceleration_velocity_twophase_explicit(double dt) noexcept { bool status = false; const double tolerance = 1.0E-15; - if (this->mass(mpm::NodePhase::nSolid) > tolerance && - this->mass(mpm::NodePhase::nLiquid) > tolerance) { + if (this->mass(mpm::NodePhase::NSolid) > tolerance && + this->mass(mpm::NodePhase::NLiquid) > tolerance) { // Compute drag force VectorDim drag_force = drag_force_coefficient_.cwiseProduct( - velocity_.col(mpm::NodePhase::nLiquid) - - velocity_.col(mpm::NodePhase::nSolid)); + velocity_.col(mpm::NodePhase::NLiquid) - + velocity_.col(mpm::NodePhase::NSolid)); // Acceleration of pore fluid (momentume balance of fluid phase) - this->acceleration_.col(mpm::NodePhase::nLiquid) = - (this->external_force_.col(mpm::NodePhase::nLiquid) + - this->internal_force_.col(mpm::NodePhase::nLiquid) - drag_force) / - this->mass_(mpm::NodePhase::nLiquid); + this->acceleration_.col(mpm::NodePhase::NLiquid) = + (this->external_force_.col(mpm::NodePhase::NLiquid) + + this->internal_force_.col(mpm::NodePhase::NLiquid) - drag_force) / + this->mass_(mpm::NodePhase::NLiquid); // Acceleration of solid skeleton (momentume balance of mixture) - this->acceleration_.col(mpm::NodePhase::nSolid) = - (this->external_force_.col(mpm::NodePhase::nMixture) + - this->internal_force_.col(mpm::NodePhase::nMixture) - - this->mass_(mpm::NodePhase::nLiquid) * - this->acceleration_.col(mpm::NodePhase::nLiquid)) / - this->mass_(mpm::NodePhase::nSolid); + this->acceleration_.col(mpm::NodePhase::NSolid) = + (this->external_force_.col(mpm::NodePhase::NMixture) + + this->internal_force_.col(mpm::NodePhase::NMixture) - + this->mass_(mpm::NodePhase::NLiquid) * + this->acceleration_.col(mpm::NodePhase::NLiquid)) / + this->mass_(mpm::NodePhase::NSolid); // Apply friction constraints this->apply_friction_constraints(dt); @@ -52,14 +52,14 @@ bool mpm::Node:: // Set a threshold for (unsigned i = 0; i < Tdim; ++i) { - if (std::abs(velocity_.col(mpm::NodePhase::nSolid)(i)) < tolerance) - velocity_.col(mpm::NodePhase::nSolid)(i) = 0.; - if (std::abs(acceleration_.col(mpm::NodePhase::nSolid)(i)) < tolerance) - acceleration_.col(mpm::NodePhase::nSolid)(i) = 0.; - if (std::abs(velocity_.col(mpm::NodePhase::nLiquid)(i)) < tolerance) - velocity_.col(mpm::NodePhase::nLiquid)(i) = 0.; - if (std::abs(acceleration_.col(mpm::NodePhase::nLiquid)(i)) < tolerance) - acceleration_.col(mpm::NodePhase::nLiquid)(i) = 0.; + if (std::abs(velocity_.col(mpm::NodePhase::NSolid)(i)) < tolerance) + velocity_.col(mpm::NodePhase::NSolid)(i) = 0.; + if (std::abs(acceleration_.col(mpm::NodePhase::NSolid)(i)) < tolerance) + acceleration_.col(mpm::NodePhase::NSolid)(i) = 0.; + if (std::abs(velocity_.col(mpm::NodePhase::NLiquid)(i)) < tolerance) + velocity_.col(mpm::NodePhase::NLiquid)(i) = 0.; + if (std::abs(acceleration_.col(mpm::NodePhase::NLiquid)(i)) < tolerance) + acceleration_.col(mpm::NodePhase::NLiquid)(i) = 0.; } status = true; } @@ -74,36 +74,36 @@ bool mpm::Node:: bool status = false; const double tolerance = 1.0E-15; - if (this->mass(mpm::NodePhase::nSolid) > tolerance && - this->mass(mpm::NodePhase::nLiquid) > tolerance) { + if (this->mass(mpm::NodePhase::NSolid) > tolerance && + this->mass(mpm::NodePhase::NLiquid) > tolerance) { // Compute drag force VectorDim drag_force = drag_force_coefficient_.cwiseProduct( - velocity_.col(mpm::NodePhase::nLiquid) - - velocity_.col(mpm::NodePhase::nSolid)); + velocity_.col(mpm::NodePhase::NLiquid) - + velocity_.col(mpm::NodePhase::NSolid)); // Unbalanced force of liquid phase auto unbalanced_force_liquid = - this->external_force_.col(mpm::NodePhase::nLiquid) + - this->internal_force_.col(mpm::NodePhase::nLiquid) - drag_force; + this->external_force_.col(mpm::NodePhase::NLiquid) + + this->internal_force_.col(mpm::NodePhase::NLiquid) - drag_force; // Acceleration of liquid phase (momentume balance of fluid phase) - this->acceleration_.col(mpm::NodePhase::nLiquid) = + this->acceleration_.col(mpm::NodePhase::NLiquid) = (unbalanced_force_liquid - damping_factor * unbalanced_force_liquid.norm() * - this->velocity_.col(mpm::NodePhase::nLiquid).cwiseSign()) / - this->mass(mpm::NodePhase::nLiquid); + this->velocity_.col(mpm::NodePhase::NLiquid).cwiseSign()) / + this->mass(mpm::NodePhase::NLiquid); // Unbalanced force of solid phase auto unbalanced_force_solid = - this->external_force_.col(mpm::NodePhase::nMixture) + - this->internal_force_.col(mpm::NodePhase::nMixture) - - this->mass_(mpm::NodePhase::nLiquid) * - this->acceleration_.col(mpm::NodePhase::nLiquid); + this->external_force_.col(mpm::NodePhase::NMixture) + + this->internal_force_.col(mpm::NodePhase::NMixture) - + this->mass_(mpm::NodePhase::NLiquid) * + this->acceleration_.col(mpm::NodePhase::NLiquid); // Acceleration of solid phase (momentume balance of mixture) - this->acceleration_.col(mpm::NodePhase::nSolid) = + this->acceleration_.col(mpm::NodePhase::NSolid) = (unbalanced_force_solid - damping_factor * unbalanced_force_solid.norm() * - this->velocity_.col(mpm::NodePhase::nSolid).cwiseSign()) / - this->mass(mpm::NodePhase::nSolid); + this->velocity_.col(mpm::NodePhase::NSolid).cwiseSign()) / + this->mass(mpm::NodePhase::NSolid); // Apply friction constraints this->apply_friction_constraints(dt); @@ -117,14 +117,14 @@ bool mpm::Node:: // Set a threshold for (unsigned i = 0; i < Tdim; ++i) { - if (std::abs(velocity_.col(mpm::NodePhase::nSolid)(i)) < tolerance) - velocity_.col(mpm::NodePhase::nSolid)(i) = 0.; - if (std::abs(acceleration_.col(mpm::NodePhase::nSolid)(i)) < tolerance) - acceleration_.col(mpm::NodePhase::nSolid)(i) = 0.; - if (std::abs(velocity_.col(mpm::NodePhase::nLiquid)(i)) < tolerance) - velocity_.col(mpm::NodePhase::nLiquid)(i) = 0.; - if (std::abs(acceleration_.col(mpm::NodePhase::nLiquid)(i)) < tolerance) - acceleration_.col(mpm::NodePhase::nLiquid)(i) = 0.; + if (std::abs(velocity_.col(mpm::NodePhase::NSolid)(i)) < tolerance) + velocity_.col(mpm::NodePhase::NSolid)(i) = 0.; + if (std::abs(acceleration_.col(mpm::NodePhase::NSolid)(i)) < tolerance) + acceleration_.col(mpm::NodePhase::NSolid)(i) = 0.; + if (std::abs(velocity_.col(mpm::NodePhase::NLiquid)(i)) < tolerance) + velocity_.col(mpm::NodePhase::NLiquid)(i) = 0.; + if (std::abs(acceleration_.col(mpm::NodePhase::NLiquid)(i)) < tolerance) + acceleration_.col(mpm::NodePhase::NLiquid)(i) = 0.; } status = true; } diff --git a/include/solvers/mpm_explicit_twophase.tcc b/include/solvers/mpm_explicit_twophase.tcc index debb780db..822e8394e 100644 --- a/include/solvers/mpm_explicit_twophase.tcc +++ b/include/solvers/mpm_explicit_twophase.tcc @@ -161,29 +161,29 @@ bool mpm::MPMExplicitTwoPhase::solve() { // MPI all reduce nodal mass for solid phase mesh_->template nodal_halo_exchange( std::bind(&mpm::NodeBase::mass, std::placeholders::_1, - mpm::NodePhase::nSolid), + mpm::NodePhase::NSolid), std::bind(&mpm::NodeBase::update_mass, std::placeholders::_1, - false, mpm::NodePhase::nSolid, std::placeholders::_2)); + false, mpm::NodePhase::NSolid, std::placeholders::_2)); // MPI all reduce nodal momentum for solid phase mesh_->template nodal_halo_exchange, Tdim>( std::bind(&mpm::NodeBase::momentum, std::placeholders::_1, - mpm::NodePhase::nSolid), + mpm::NodePhase::NSolid), std::bind(&mpm::NodeBase::update_momentum, - std::placeholders::_1, false, mpm::NodePhase::nSolid, + std::placeholders::_1, false, mpm::NodePhase::NSolid, std::placeholders::_2)); // MPI all reduce nodal mass for liquid phase mesh_->template nodal_halo_exchange( std::bind(&mpm::NodeBase::mass, std::placeholders::_1, - mpm::NodePhase::nLiquid), + mpm::NodePhase::NLiquid), std::bind(&mpm::NodeBase::update_mass, std::placeholders::_1, - false, mpm::NodePhase::nLiquid, std::placeholders::_2)); + false, mpm::NodePhase::NLiquid, std::placeholders::_2)); // MPI all reduce nodal momentum for liquid phase mesh_->template nodal_halo_exchange, Tdim>( std::bind(&mpm::NodeBase::momentum, std::placeholders::_1, - mpm::NodePhase::nLiquid), + mpm::NodePhase::NLiquid), std::bind(&mpm::NodeBase::update_momentum, - std::placeholders::_1, false, mpm::NodePhase::nLiquid, + std::placeholders::_1, false, mpm::NodePhase::NLiquid, std::placeholders::_2)); } #endif @@ -263,31 +263,31 @@ bool mpm::MPMExplicitTwoPhase::solve() { // MPI all reduce external force of mixture mesh_->template nodal_halo_exchange, Tdim>( std::bind(&mpm::NodeBase::external_force, std::placeholders::_1, - mpm::NodePhase::nMixture), + mpm::NodePhase::NMixture), std::bind(&mpm::NodeBase::update_external_force, - std::placeholders::_1, false, mpm::NodePhase::nMixture, + std::placeholders::_1, false, mpm::NodePhase::NMixture, std::placeholders::_2)); // MPI all reduce external force of pore fluid mesh_->template nodal_halo_exchange, Tdim>( std::bind(&mpm::NodeBase::external_force, std::placeholders::_1, - mpm::NodePhase::nLiquid), + mpm::NodePhase::NLiquid), std::bind(&mpm::NodeBase::update_external_force, - std::placeholders::_1, false, mpm::NodePhase::nLiquid, + std::placeholders::_1, false, mpm::NodePhase::NLiquid, std::placeholders::_2)); // MPI all reduce internal force of mixture mesh_->template nodal_halo_exchange, Tdim>( std::bind(&mpm::NodeBase::internal_force, std::placeholders::_1, - mpm::NodePhase::nMixture), + mpm::NodePhase::NMixture), std::bind(&mpm::NodeBase::update_internal_force, - std::placeholders::_1, false, mpm::NodePhase::nMixture, + std::placeholders::_1, false, mpm::NodePhase::NMixture, std::placeholders::_2)); // MPI all reduce internal force of pore liquid mesh_->template nodal_halo_exchange, Tdim>( std::bind(&mpm::NodeBase::internal_force, std::placeholders::_1, - mpm::NodePhase::nLiquid), + mpm::NodePhase::NLiquid), std::bind(&mpm::NodeBase::update_internal_force, - std::placeholders::_1, false, mpm::NodePhase::nLiquid, + std::placeholders::_1, false, mpm::NodePhase::NLiquid, std::placeholders::_2)); // MPI all reduce drag force diff --git a/tests/nodes/node_test.cc b/tests/nodes/node_test.cc index 7552300eb..8b2d4968d 100644 --- a/tests/nodes/node_test.cc +++ b/tests/nodes/node_test.cc @@ -159,13 +159,13 @@ TEST_CASE("Node is checked for 1D case", "[node][1D]") { // Check pressure constraints SECTION("Check nodal pressure constraints") { // Check assign pressure constraint - REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::nSolid, 8000, + REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::NSolid, 8000, nullptr) == true); // Check apply pressure constraint REQUIRE_NOTHROW( - node->apply_pressure_constraint(mpm::NodePhase::nSolid)); + node->apply_pressure_constraint(mpm::NodePhase::NSolid)); // Check pressure - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(8000).epsilon(Tolerance)); } } @@ -636,13 +636,13 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { // Check pressure constraints SECTION("Check nodal pressure constraints") { // Check assign pressure constraint - REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::nSolid, 8000, + REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::NSolid, 8000, nullptr) == true); // Check apply pressure constraint REQUIRE_NOTHROW( - node->apply_pressure_constraint(mpm::NodePhase::nSolid)); + node->apply_pressure_constraint(mpm::NodePhase::NSolid)); // Check pressure - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(8000).epsilon(Tolerance)); } } @@ -1250,13 +1250,13 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { // Check pressure constraints SECTION("Check nodal pressure constraints") { // Check assign pressure constraint - REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::nSolid, 8000, + REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::NSolid, 8000, nullptr) == true); // Check apply pressure constraint REQUIRE_NOTHROW( - node->apply_pressure_constraint(mpm::NodePhase::nSolid)); + node->apply_pressure_constraint(mpm::NodePhase::NSolid)); // Check pressure - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(8000).epsilon(Tolerance)); } } diff --git a/tests/nodes/node_twophase_test.cc b/tests/nodes/node_twophase_test.cc index eebb0182c..d5fa22f1f 100644 --- a/tests/nodes/node_twophase_test.cc +++ b/tests/nodes/node_twophase_test.cc @@ -121,111 +121,111 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { std::make_shared>(id, coords); // Check mass - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); double solid_mass = 100.5; double liquid_mass = 200.5; // Update mass to 100.5 and 200.5 REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(true, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(true, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(100.5).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(200.5).epsilon(Tolerance)); // Update mass to 201 and 401 REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(true, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(true, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(201.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(401.0).epsilon(Tolerance)); // Assign mass to 100 and 200 solid_mass = 100.; liquid_mass = 200.; REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(false, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(false, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(100.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(200.0).epsilon(Tolerance)); SECTION("Check nodal pressure") { // Check pressure - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); double pressure = 1000.7; double pore_pressure = 2000.7; // Update pressure to 1000.7 and 2000.7 - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nSolid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NSolid, solid_mass * pressure)); - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nLiquid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NLiquid, liquid_mass * pore_pressure)); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(1000.7).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(2000.7).epsilon(Tolerance)); // Update pressure to 2001.4 and 4001.4 - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nSolid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NSolid, solid_mass * pressure)); - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nLiquid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NLiquid, liquid_mass * pore_pressure)); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(2001.4).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(4001.4).epsilon(Tolerance)); // Assign pressure to 1000 and 2000 pressure = 1000.; pore_pressure = 2000.; - node->assign_pressure(mpm::NodePhase::nSolid, pressure); - node->assign_pressure(mpm::NodePhase::nLiquid, pore_pressure); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + node->assign_pressure(mpm::NodePhase::NSolid, pressure); + node->assign_pressure(mpm::NodePhase::NLiquid, pore_pressure); + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(1000.0).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(2000.0).epsilon(Tolerance)); // Assign mass to 0 solid_mass = 0.; liquid_mass = 0.; REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(false, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(false, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); // Try to update pressure to 2000, should throw and keep to 1000. - node->assign_pressure(mpm::NodePhase::nSolid, pressure); - node->assign_pressure(mpm::NodePhase::nLiquid, pore_pressure); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + node->assign_pressure(mpm::NodePhase::NSolid, pressure); + node->assign_pressure(mpm::NodePhase::NLiquid, pore_pressure); + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(1000.0).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(2000.0).epsilon(Tolerance)); // Check pressure constraints SECTION("Check nodal pressure constraints") { // Check assign pressure constraint - REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::nSolid, 8000, + REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::NSolid, 8000, nullptr) == true); - REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::nLiquid, 7000, + REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::NLiquid, 7000, nullptr) == true); // Check apply pressure constraint REQUIRE_NOTHROW( - node->apply_pressure_constraint(mpm::NodePhase::nSolid)); + node->apply_pressure_constraint(mpm::NodePhase::NSolid)); REQUIRE_NOTHROW( - node->apply_pressure_constraint(mpm::NodePhase::nLiquid)); + node->apply_pressure_constraint(mpm::NodePhase::NLiquid)); // Check pressure - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(8000).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(7000).epsilon(Tolerance)); } } @@ -237,45 +237,45 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { // Check current external force is zero for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Update force to 10.0 REQUIRE_NOTHROW( - node->update_external_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::nLiquid, + node->update_external_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } // Update force to 20.0 REQUIRE_NOTHROW( - node->update_external_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::nLiquid, + node->update_external_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(20.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Assign force as 10.0 REQUIRE_NOTHROW( - node->update_external_force(false, mpm::NodePhase::nMixture, force)); + node->update_external_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_external_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } @@ -283,7 +283,7 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { // Set external force to zero force.setZero(); REQUIRE_NOTHROW(node->update_external_force( - false, mpm::NodePhase::nMixture, force)); + false, mpm::NodePhase::NMixture, force)); // concentrated force std::shared_ptr ffunction = nullptr; @@ -291,38 +291,38 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { const unsigned Direction = 0; // Check external force for (unsigned i = 0; i < Dim; ++i) - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->assign_concentrated_force(mpm::NodePhase::nMixture, + REQUIRE(node->assign_concentrated_force(mpm::NodePhase::NMixture, Direction, concentrated_force, ffunction) == true); double current_time = 0.0; - node->apply_concentrated_force(mpm::NodePhase::nMixture, current_time); + node->apply_concentrated_force(mpm::NodePhase::NMixture, current_time); for (unsigned i = 0; i < Dim; ++i) { if (i == Direction) - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(concentrated_force).epsilon(Tolerance)); else - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); } // Check for incorrect direction / phase const unsigned wrong_dir = 4; - REQUIRE(node->assign_concentrated_force(mpm::NodePhase::nMixture, + REQUIRE(node->assign_concentrated_force(mpm::NodePhase::NMixture, wrong_dir, concentrated_force, ffunction) == false); // Check again to ensure value hasn't been updated for (unsigned i = 0; i < Dim; ++i) { if (i == Direction) - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(concentrated_force).epsilon(Tolerance)); else - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); } } @@ -335,45 +335,45 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { // Check current internal force is zero for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Update force to 10.0 REQUIRE_NOTHROW( - node->update_internal_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::nLiquid, + node->update_internal_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } // Update force to 20.0 REQUIRE_NOTHROW( - node->update_internal_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::nLiquid, + node->update_internal_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(20.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Assign force as 10.0 REQUIRE_NOTHROW( - node->update_internal_force(false, mpm::NodePhase::nMixture, force)); + node->update_internal_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_internal_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } } @@ -422,12 +422,12 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { double liquid_mass = 100.; // Update mass to 100.5 REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nSolid, solid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(false, mpm::NodePhase::NSolid, solid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(solid_mass).epsilon(Tolerance)); REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + node->update_mass(false, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(liquid_mass).epsilon(Tolerance)); // Check internal force @@ -436,14 +436,14 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { for (unsigned i = 0; i < force.size(); ++i) force(i) = 10. * i; // Update force to 10.0 REQUIRE_NOTHROW( - node->update_internal_force(false, mpm::NodePhase::nMixture, force)); + node->update_internal_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_internal_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); // Internal force for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(force(i)).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(0.5 * force(i)).epsilon(Tolerance)); } @@ -451,13 +451,13 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { for (unsigned i = 0; i < force.size(); ++i) force(i) = 5. * i; // Update force to 10.0 REQUIRE_NOTHROW( - node->update_external_force(false, mpm::NodePhase::nMixture, force)); + node->update_external_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_external_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(force(i)).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(0.5 * force(i)).epsilon(Tolerance)); } @@ -482,9 +482,9 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { solid_acceleration << 0.; for (unsigned i = 0; i < solid_acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(solid_acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(liquid_acceleration(i)).epsilon(Tolerance)); } @@ -492,9 +492,9 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { Eigen::Matrix solid_velocity = solid_acceleration * dt; Eigen::Matrix liquid_velocity = liquid_acceleration * dt; for (unsigned i = 0; i < solid_velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(solid_velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(liquid_velocity(i)).epsilon(Tolerance)); } @@ -508,9 +508,9 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { solid_acceleration[0] = 0.5 * solid_acceleration[0]; liquid_acceleration[0] = 0.5 * liquid_acceleration[0]; for (unsigned i = 0; i < solid_acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(solid_acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(liquid_acceleration(i)).epsilon(Tolerance)); } @@ -522,9 +522,9 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { solid_acceleration[0] = 0.; liquid_acceleration[0] = 0.; for (unsigned i = 0; i < solid_acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(solid_acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(liquid_acceleration(i)).epsilon(Tolerance)); } @@ -538,9 +538,9 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { solid_velocity[0] = 10.5; liquid_velocity[0] = 10.5; for (unsigned i = 0; i < solid_velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(solid_velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(liquid_velocity(i)).epsilon(Tolerance)); } @@ -548,19 +548,19 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { solid_acceleration[0] = 0.; liquid_acceleration[0] = 0.; for (unsigned i = 0; i < solid_acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(solid_acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(liquid_acceleration(i)).epsilon(Tolerance)); } // Exception check when mass is zero // Update mass to 0. - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nSolid, 0.)); - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nLiquid, 0.)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NSolid, 0.)); + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NLiquid, 0.)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.).epsilon(Tolerance)); REQUIRE(node->compute_acceleration_velocity_twophase_explicit(dt) == false); @@ -573,82 +573,82 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { // Check initial momentum for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Check update momentum to 10 REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nSolid, momentum)); + node->update_momentum(true, mpm::NodePhase::NSolid, momentum)); REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nLiquid, momentum)); + node->update_momentum(true, mpm::NodePhase::NLiquid, momentum)); for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Check update momentum to 20 REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nSolid, momentum)); + node->update_momentum(true, mpm::NodePhase::NSolid, momentum)); REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nLiquid, momentum)); + node->update_momentum(true, mpm::NodePhase::NLiquid, momentum)); for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(20.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(20.).epsilon(Tolerance)); } // Check assign momentum to 10 REQUIRE_NOTHROW( - node->update_momentum(false, mpm::NodePhase::nSolid, momentum)); + node->update_momentum(false, mpm::NodePhase::NSolid, momentum)); REQUIRE_NOTHROW( - node->update_momentum(false, mpm::NodePhase::nLiquid, momentum)); + node->update_momentum(false, mpm::NodePhase::NLiquid, momentum)); for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Check zero velocity for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Check mass double mass = 0.; - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nSolid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NSolid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nLiquid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NLiquid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); // Compute and check velocity this should throw zero mass node->compute_velocity(); mass = 100.; // Update mass to 100. - REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::nSolid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::NSolid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(100.).epsilon(Tolerance)); - REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::nLiquid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::NLiquid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(100.).epsilon(Tolerance)); // Compute and check velocity node->compute_velocity(); for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(0.1).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(0.1).epsilon(Tolerance)); } @@ -661,9 +661,9 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { Eigen::Matrix velocity; velocity << 0.1; for (unsigned i = 0; i < velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(velocity(i)).epsilon(Tolerance)); } @@ -673,9 +673,9 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { // Check apply constraints velocity << 10.5; for (unsigned i = 0; i < velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(velocity(i)).epsilon(Tolerance)); } } @@ -686,21 +686,21 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { for (unsigned i = 0; i < acceleration.size(); ++i) acceleration(i) = 5.; for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } - REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::nSolid, + REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::NSolid, acceleration)); - REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::nLiquid, + REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::NLiquid, 0.5 * acceleration)); for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(5.).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(2.5).epsilon(Tolerance)); } @@ -713,9 +713,9 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { acceleration.resize(Dim); acceleration << 5.; for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(0.5 * acceleration(i)).epsilon(Tolerance)); } @@ -725,9 +725,9 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") { // Check apply constraints acceleration << 0.0; for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } } @@ -866,135 +866,135 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { std::make_shared>(id, coords); // Check mass - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); double solid_mass = 100.5; double liquid_mass = 200.5; // Update mass to 100.5 and 200.5 REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(true, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(true, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(100.5).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(200.5).epsilon(Tolerance)); // Update mass to 201 and 401 REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(true, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(true, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(201.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(401.0).epsilon(Tolerance)); // Assign mass to 100 and 200 solid_mass = 100.; liquid_mass = 200.; REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(false, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(false, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(100.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(200.0).epsilon(Tolerance)); SECTION("Check nodal pressure") { // Check pressure - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); double pressure = 1000.7; double pore_pressure = 2000.7; // Update pressure to 1000.7 and 2000.7 - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nSolid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NSolid, solid_mass * pressure)); - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nLiquid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NLiquid, liquid_mass * pore_pressure)); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(1000.7).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(2000.7).epsilon(Tolerance)); // Update pressure to 2001.4 and 4001.4 - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nSolid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NSolid, solid_mass * pressure)); - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nLiquid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NLiquid, liquid_mass * pore_pressure)); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(2001.4).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(4001.4).epsilon(Tolerance)); // Assign pressure to 1000 and 2000 pressure = 1000.; pore_pressure = 2000.; - node->assign_pressure(mpm::NodePhase::nSolid, pressure); - node->assign_pressure(mpm::NodePhase::nLiquid, pore_pressure); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + node->assign_pressure(mpm::NodePhase::NSolid, pressure); + node->assign_pressure(mpm::NodePhase::NLiquid, pore_pressure); + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(1000.0).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(2000.0).epsilon(Tolerance)); // Assign mass to 0 solid_mass = 0.; liquid_mass = 0.; REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(false, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(false, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); // Try to update pressure to 2000, should throw and keep to 1000. - node->assign_pressure(mpm::NodePhase::nSolid, pressure); - node->assign_pressure(mpm::NodePhase::nLiquid, pore_pressure); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + node->assign_pressure(mpm::NodePhase::NSolid, pressure); + node->assign_pressure(mpm::NodePhase::NLiquid, pore_pressure); + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(1000.0).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(2000.0).epsilon(Tolerance)); // Check pressure constraints SECTION("Check nodal pressure constraints") { // Check assign pressure constraint - REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::nSolid, 8000, + REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::NSolid, 8000, nullptr) == true); - REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::nLiquid, 7000, + REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::NLiquid, 7000, nullptr) == true); // Check apply pressure constraint REQUIRE_NOTHROW( - node->apply_pressure_constraint(mpm::NodePhase::nSolid)); + node->apply_pressure_constraint(mpm::NodePhase::NSolid)); REQUIRE_NOTHROW( - node->apply_pressure_constraint(mpm::NodePhase::nLiquid)); + node->apply_pressure_constraint(mpm::NodePhase::NLiquid)); // Check pressure - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(8000).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(7000).epsilon(Tolerance)); } } SECTION("Check volume") { // Check volume - REQUIRE(node->volume(mpm::NodePhase::nMixture) == + REQUIRE(node->volume(mpm::NodePhase::NMixture) == Approx(0.0).epsilon(Tolerance)); double volume = 100.5; // Update volume to 100.5 REQUIRE_NOTHROW( - node->update_volume(true, mpm::NodePhase::nMixture, volume)); - REQUIRE(node->volume(mpm::NodePhase::nMixture) == + node->update_volume(true, mpm::NodePhase::NMixture, volume)); + REQUIRE(node->volume(mpm::NodePhase::NMixture) == Approx(100.5).epsilon(Tolerance)); // Update volume to 201 REQUIRE_NOTHROW( - node->update_volume(true, mpm::NodePhase::nMixture, volume)); - REQUIRE(node->volume(mpm::NodePhase::nMixture) == + node->update_volume(true, mpm::NodePhase::NMixture, volume)); + REQUIRE(node->volume(mpm::NodePhase::NMixture) == Approx(201.0).epsilon(Tolerance)); // Assign volume to 100 volume = 100.; REQUIRE_NOTHROW( - node->update_volume(false, mpm::NodePhase::nMixture, volume)); - REQUIRE(node->volume(mpm::NodePhase::nMixture) == + node->update_volume(false, mpm::NodePhase::NMixture, volume)); + REQUIRE(node->volume(mpm::NodePhase::NMixture) == Approx(100.0).epsilon(Tolerance)); } @@ -1005,45 +1005,45 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { // Check current external force is zero for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Update force to 10.0 REQUIRE_NOTHROW( - node->update_external_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::nLiquid, + node->update_external_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } // Update force to 20.0 REQUIRE_NOTHROW( - node->update_external_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::nLiquid, + node->update_external_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(20.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Assign force as 10.0 REQUIRE_NOTHROW( - node->update_external_force(false, mpm::NodePhase::nMixture, force)); + node->update_external_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_external_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } @@ -1051,7 +1051,7 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { // Set external force to zero force.setZero(); REQUIRE_NOTHROW(node->update_external_force( - false, mpm::NodePhase::nMixture, force)); + false, mpm::NodePhase::NMixture, force)); // Concentrated force std::shared_ptr ffunction = nullptr; @@ -1059,38 +1059,38 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { const unsigned Direction = 0; // Check traction for (unsigned i = 0; i < Dim; ++i) - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->assign_concentrated_force(mpm::NodePhase::nMixture, + REQUIRE(node->assign_concentrated_force(mpm::NodePhase::NMixture, Direction, concentrated_force, ffunction) == true); double current_time = 0.0; - node->apply_concentrated_force(mpm::NodePhase::nMixture, current_time); + node->apply_concentrated_force(mpm::NodePhase::NMixture, current_time); for (unsigned i = 0; i < Dim; ++i) { if (i == Direction) - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(concentrated_force).epsilon(Tolerance)); else - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); } // Check for incorrect direction / phase const unsigned wrong_dir = 4; - REQUIRE(node->assign_concentrated_force(mpm::NodePhase::nMixture, + REQUIRE(node->assign_concentrated_force(mpm::NodePhase::NMixture, wrong_dir, concentrated_force, ffunction) == false); // Check again to ensure value hasn't been updated for (unsigned i = 0; i < Dim; ++i) { if (i == Direction) - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(concentrated_force).epsilon(Tolerance)); else - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); } } @@ -1103,45 +1103,45 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { // Check current internal force is zero for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Update force to 10.0 REQUIRE_NOTHROW( - node->update_internal_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::nLiquid, + node->update_internal_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } // Update force to 20.0 REQUIRE_NOTHROW( - node->update_internal_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::nLiquid, + node->update_internal_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(20.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Assign force as 10.0 REQUIRE_NOTHROW( - node->update_internal_force(false, mpm::NodePhase::nMixture, force)); + node->update_internal_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_internal_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } } @@ -1190,12 +1190,12 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { double liquid_mass = 100.; // Update mass to 100. REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(false, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(false, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(solid_mass).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(liquid_mass).epsilon(Tolerance)); // Check internal force @@ -1204,14 +1204,14 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { for (unsigned i = 0; i < force.size(); ++i) force(i) = 10. * i; // Update force to 10.0 REQUIRE_NOTHROW( - node->update_internal_force(false, mpm::NodePhase::nMixture, force)); + node->update_internal_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_internal_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); // Internal force for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(force(i)).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(0.5 * force(i)).epsilon(Tolerance)); } @@ -1219,13 +1219,13 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { for (unsigned i = 0; i < force.size(); ++i) force(i) = 5. * i; // Update force to 10.0 REQUIRE_NOTHROW( - node->update_external_force(false, mpm::NodePhase::nMixture, force)); + node->update_external_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_external_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(force(i)).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(0.5 * force(i)).epsilon(Tolerance)); } @@ -1250,9 +1250,9 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { solid_acceleration << 0., 0.075; for (unsigned i = 0; i < solid_acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(solid_acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(liquid_acceleration(i)).epsilon(Tolerance)); } @@ -1260,9 +1260,9 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { Eigen::Matrix solid_velocity = solid_acceleration * dt; Eigen::Matrix liquid_velocity = liquid_acceleration * dt; for (unsigned i = 0; i < solid_velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(solid_velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(liquid_velocity(i)).epsilon(Tolerance)); } @@ -1278,9 +1278,9 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { solid_velocity << 10.5, 0.03; liquid_velocity << 20.5, 1.03; for (unsigned i = 0; i < solid_velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(solid_velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(liquid_velocity(i)).epsilon(Tolerance)); } @@ -1288,9 +1288,9 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { solid_acceleration.setZero(); liquid_acceleration.setZero(); for (unsigned i = 0; i < solid_acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(solid_acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(liquid_acceleration(i)).epsilon(Tolerance)); } @@ -1302,18 +1302,18 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { solid_acceleration << 0., 0.; liquid_acceleration << 0., 0.; for (unsigned i = 0; i < solid_acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(solid_acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(liquid_acceleration(i)).epsilon(Tolerance)); } // Exception check when mass is zero - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nSolid, 0.)); - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nLiquid, 0.)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NSolid, 0.)); + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NLiquid, 0.)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.).epsilon(Tolerance)); REQUIRE(node->compute_acceleration_velocity_twophase_explicit(dt) == false); @@ -1329,82 +1329,82 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { // Check initial momentum for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Check update momentum to 10 REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nSolid, momentum)); + node->update_momentum(true, mpm::NodePhase::NSolid, momentum)); REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nLiquid, momentum)); + node->update_momentum(true, mpm::NodePhase::NLiquid, momentum)); for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Check update momentum to 20 REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nSolid, momentum)); + node->update_momentum(true, mpm::NodePhase::NSolid, momentum)); REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nLiquid, momentum)); + node->update_momentum(true, mpm::NodePhase::NLiquid, momentum)); for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(20.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(20.).epsilon(Tolerance)); } // Check assign momentum to 10 REQUIRE_NOTHROW( - node->update_momentum(false, mpm::NodePhase::nSolid, momentum)); + node->update_momentum(false, mpm::NodePhase::NSolid, momentum)); REQUIRE_NOTHROW( - node->update_momentum(false, mpm::NodePhase::nLiquid, momentum)); + node->update_momentum(false, mpm::NodePhase::NLiquid, momentum)); for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Check zero velocity for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Check mass double mass = 0.; - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nSolid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NSolid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nLiquid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NLiquid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); // Compute and check velocity this should throw zero mass node->compute_velocity(); mass = 100.; // Update mass to 100.5 - REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::nSolid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::NSolid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(100.).epsilon(Tolerance)); - REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::nLiquid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::NLiquid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(100.).epsilon(Tolerance)); // Compute and check velocity node->compute_velocity(); for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(0.1).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(0.1).epsilon(Tolerance)); } @@ -1413,20 +1413,20 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { for (unsigned i = 0; i < acceleration.size(); ++i) acceleration(i) = 5.; for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } - REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::nSolid, + REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::NSolid, acceleration)); - REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::nLiquid, + REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::NLiquid, acceleration)); for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(5.).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } @@ -1439,18 +1439,18 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { Eigen::Matrix velocity; velocity << 0.1, 0.1; for (unsigned i = 0; i < velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(velocity(i)).epsilon(Tolerance)); } // Check acceleration before constraints acceleration << 5., 5.; for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } @@ -1466,17 +1466,17 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { // Check apply constraints velocity << -12.5, 0.1; for (unsigned i = 0; i < velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(velocity(i)).epsilon(Tolerance)); } acceleration << 0., 5.; for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } } @@ -1501,35 +1501,35 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { // Check applied velocity constraints in the global coordinates velocity << -9.583478335521184, -8.025403099849004; for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(velocity(i)).epsilon(Tolerance)); } // Check that the velocity is as specified in local coordinate REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nSolid))(0) == + node->velocity(mpm::NodePhase::NSolid))(0) == Approx(-12.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nLiquid))(0) == + node->velocity(mpm::NodePhase::NLiquid))(0) == Approx(-12.5).epsilon(Tolerance)); // Check applied constraints on acceleration in the global coordinates acceleration << -0.396139826697847, 0.472101061636807; for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } // Check that the acceleration is 0 in local coordinate REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nSolid))(0) == + node->acceleration(mpm::NodePhase::NSolid))(0) == Approx(0).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nLiquid))(0) == + node->acceleration(mpm::NodePhase::NLiquid))(0) == Approx(0).epsilon(Tolerance)); } @@ -1555,47 +1555,47 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { // Check applied velocity constraints in the global coordinates velocity << -14.311308834766370, 2.772442864323454; for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(velocity(i)).epsilon(Tolerance)); } // Check that the velocity is as specified in local coordinate REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nSolid))(0) == + node->velocity(mpm::NodePhase::NSolid))(0) == Approx(-12.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nSolid))(1) == + node->velocity(mpm::NodePhase::NSolid))(1) == Approx(7.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nLiquid))(0) == + node->velocity(mpm::NodePhase::NLiquid))(0) == Approx(-12.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nLiquid))(1) == + node->velocity(mpm::NodePhase::NLiquid))(1) == Approx(7.5).epsilon(Tolerance)); // Check applied constraints on acceleration in the global coordinates acceleration << 0, 0; for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } // Check that the acceleration is 0 in local coordinate REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nSolid))(0) == + node->acceleration(mpm::NodePhase::NSolid))(0) == Approx(0).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nSolid))(1) == + node->acceleration(mpm::NodePhase::NSolid))(1) == Approx(0).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nLiquid))(0) == + node->acceleration(mpm::NodePhase::NLiquid))(0) == Approx(0).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nLiquid))(1) == + node->acceleration(mpm::NodePhase::NLiquid))(1) == Approx(0).epsilon(Tolerance)); } @@ -1611,7 +1611,7 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { // Check apply constraints acceleration << 4., 5.; for (unsigned i = 0; i < acceleration.size(); ++i) - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } @@ -1633,14 +1633,14 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") { // Check applied constraints on acceleration in the global coordinates acceleration << 4.905579787672637, 4.920772034660430; for (unsigned i = 0; i < Dim; ++i) - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); // Check the acceleration in local coordinate acceleration << 6.920903430595146, 0.616284167162194; for (unsigned i = 0; i < Dim; ++i) REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nSolid))(i) == + node->acceleration(mpm::NodePhase::NSolid))(i) == Approx(acceleration(i)).epsilon(Tolerance)); } } @@ -1779,111 +1779,111 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { std::make_shared>(id, coords); // Check mass - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); double solid_mass = 100.5; double liquid_mass = 200.5; // Update mass to 100.5 and 200.5 REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(true, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(true, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(100.5).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(200.5).epsilon(Tolerance)); // Update mass to 201 and 401 REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(true, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(true, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(true, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(201.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(401.0).epsilon(Tolerance)); // Assign mass to 100 and 200 solid_mass = 100.; liquid_mass = 200.; REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(false, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(false, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(100.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(200.0).epsilon(Tolerance)); SECTION("Check nodal pressure") { // Check pressure - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); double pressure = 1000.7; double pore_pressure = 2000.7; // Update pressure to 1000.7 and 2000.7 - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nSolid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NSolid, solid_mass * pressure)); - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nLiquid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NLiquid, liquid_mass * pore_pressure)); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(1000.7).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(2000.7).epsilon(Tolerance)); // Update pressure to 2001.4 and 4001.4 - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nSolid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NSolid, solid_mass * pressure)); - REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::nLiquid, + REQUIRE_NOTHROW(node->update_mass_pressure(mpm::NodePhase::NLiquid, liquid_mass * pore_pressure)); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(2001.4).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(4001.4).epsilon(Tolerance)); // Assign pressure to 1000 and 2000 pressure = 1000.; pore_pressure = 2000.; - node->assign_pressure(mpm::NodePhase::nSolid, pressure); - node->assign_pressure(mpm::NodePhase::nLiquid, pore_pressure); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + node->assign_pressure(mpm::NodePhase::NSolid, pressure); + node->assign_pressure(mpm::NodePhase::NLiquid, pore_pressure); + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(1000.0).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(2000.0).epsilon(Tolerance)); // Assign mass to 0 solid_mass = 0.; liquid_mass = 0.; REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(false, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(false, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); // Try to update pressure to 2000, should throw and keep to 1000. - node->assign_pressure(mpm::NodePhase::nSolid, pressure); - node->assign_pressure(mpm::NodePhase::nLiquid, pore_pressure); - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + node->assign_pressure(mpm::NodePhase::NSolid, pressure); + node->assign_pressure(mpm::NodePhase::NLiquid, pore_pressure); + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(1000.0).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(2000.0).epsilon(Tolerance)); // Check pressure constraints SECTION("Check nodal pressure constraints") { // Check assign pressure constraint - REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::nSolid, 8000, + REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::NSolid, 8000, nullptr) == true); - REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::nLiquid, 7000, + REQUIRE(node->assign_pressure_constraint(mpm::NodePhase::NLiquid, 7000, nullptr) == true); // Check apply pressure constraint REQUIRE_NOTHROW( - node->apply_pressure_constraint(mpm::NodePhase::nSolid)); + node->apply_pressure_constraint(mpm::NodePhase::NSolid)); REQUIRE_NOTHROW( - node->apply_pressure_constraint(mpm::NodePhase::nLiquid)); + node->apply_pressure_constraint(mpm::NodePhase::NLiquid)); // Check pressure - REQUIRE(node->pressure(mpm::NodePhase::nSolid) == + REQUIRE(node->pressure(mpm::NodePhase::NSolid) == Approx(8000).epsilon(Tolerance)); - REQUIRE(node->pressure(mpm::NodePhase::nLiquid) == + REQUIRE(node->pressure(mpm::NodePhase::NLiquid) == Approx(7000).epsilon(Tolerance)); } } @@ -1895,45 +1895,45 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { // Check current external force is zero for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Update force to 10.0 REQUIRE_NOTHROW( - node->update_external_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::nLiquid, + node->update_external_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } // Update force to 20.0 REQUIRE_NOTHROW( - node->update_external_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::nLiquid, + node->update_external_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_external_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(20.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Assign force as 10.0 REQUIRE_NOTHROW( - node->update_external_force(false, mpm::NodePhase::nMixture, force)); + node->update_external_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_external_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } @@ -1941,7 +1941,7 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { // Set external force to zero force.setZero(); REQUIRE_NOTHROW(node->update_external_force( - false, mpm::NodePhase::nMixture, force)); + false, mpm::NodePhase::NMixture, force)); // Concentrated force std::shared_ptr ffunction = nullptr; @@ -1949,38 +1949,38 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { const unsigned Direction = 0; // Check traction for (unsigned i = 0; i < Dim; ++i) - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->assign_concentrated_force(mpm::NodePhase::nMixture, + REQUIRE(node->assign_concentrated_force(mpm::NodePhase::NMixture, Direction, concentrated_force, ffunction) == true); double current_time = 0.0; - node->apply_concentrated_force(mpm::NodePhase::nMixture, current_time); + node->apply_concentrated_force(mpm::NodePhase::NMixture, current_time); for (unsigned i = 0; i < Dim; ++i) { if (i == Direction) - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(concentrated_force).epsilon(Tolerance)); else - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); } // Check for incorrect direction / phase const unsigned wrong_dir = 4; - REQUIRE(node->assign_concentrated_force(mpm::NodePhase::nMixture, + REQUIRE(node->assign_concentrated_force(mpm::NodePhase::NMixture, wrong_dir, concentrated_force, ffunction) == false); // Check again to ensure value hasn't been updated for (unsigned i = 0; i < Dim; ++i) { if (i == Direction) - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(concentrated_force).epsilon(Tolerance)); else - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); } } @@ -1993,45 +1993,45 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { // Check current internal force is zero for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Update force to 10.0 REQUIRE_NOTHROW( - node->update_internal_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::nLiquid, + node->update_internal_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } // Update force to 20.0 REQUIRE_NOTHROW( - node->update_internal_force(true, mpm::NodePhase::nMixture, force)); - REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::nLiquid, + node->update_internal_force(true, mpm::NodePhase::NMixture, force)); + REQUIRE_NOTHROW(node->update_internal_force(true, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(20.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Assign force as 10.0 REQUIRE_NOTHROW( - node->update_internal_force(false, mpm::NodePhase::nMixture, force)); + node->update_internal_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_internal_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } } @@ -2080,12 +2080,12 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { double liquid_mass = 100.; // Update mass to 100. REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nSolid, solid_mass)); + node->update_mass(false, mpm::NodePhase::NSolid, solid_mass)); REQUIRE_NOTHROW( - node->update_mass(false, mpm::NodePhase::nLiquid, liquid_mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + node->update_mass(false, mpm::NodePhase::NLiquid, liquid_mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(solid_mass).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(liquid_mass).epsilon(Tolerance)); // Check internal force @@ -2094,14 +2094,14 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { for (unsigned i = 0; i < force.size(); ++i) force(i) = 10. * i; // Update force to 10.0 REQUIRE_NOTHROW( - node->update_internal_force(false, mpm::NodePhase::nMixture, force)); + node->update_internal_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_internal_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); // Internal force for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->internal_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NMixture)(i) == Approx(force(i)).epsilon(Tolerance)); - REQUIRE(node->internal_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->internal_force(mpm::NodePhase::NLiquid)(i) == Approx(0.5 * force(i)).epsilon(Tolerance)); } @@ -2109,13 +2109,13 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { for (unsigned i = 0; i < force.size(); ++i) force(i) = 5. * i; // Update force to 10.0 REQUIRE_NOTHROW( - node->update_external_force(false, mpm::NodePhase::nMixture, force)); + node->update_external_force(false, mpm::NodePhase::NMixture, force)); REQUIRE_NOTHROW(node->update_external_force( - false, mpm::NodePhase::nLiquid, 0.5 * force)); + false, mpm::NodePhase::NLiquid, 0.5 * force)); for (unsigned i = 0; i < force.size(); ++i) { - REQUIRE(node->external_force(mpm::NodePhase::nMixture)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NMixture)(i) == Approx(force(i)).epsilon(Tolerance)); - REQUIRE(node->external_force(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->external_force(mpm::NodePhase::NLiquid)(i) == Approx(0.5 * force(i)).epsilon(Tolerance)); } @@ -2141,9 +2141,9 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { // Check acceleration for (unsigned i = 0; i < solid_acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(solid_acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(liquid_acceleration(i)).epsilon(Tolerance)); } @@ -2151,9 +2151,9 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { Eigen::Matrix solid_velocity = solid_acceleration * dt; Eigen::Matrix liquid_velocity = liquid_acceleration * dt; for (unsigned i = 0; i < solid_velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(solid_velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(liquid_velocity(i)).epsilon(Tolerance)); } @@ -2171,9 +2171,9 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { solid_velocity << 10.5, 0.03, 5.13; liquid_velocity << 20.5, 1.03, 7.13; for (unsigned i = 0; i < solid_velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(solid_velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(liquid_velocity(i)).epsilon(Tolerance)); } @@ -2181,9 +2181,9 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { solid_acceleration.setZero(); liquid_acceleration.setZero(); for (unsigned i = 0; i < solid_acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(solid_acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(liquid_acceleration(i)).epsilon(Tolerance)); } @@ -2195,9 +2195,9 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { solid_acceleration << 0., 0., 0.; liquid_acceleration << 0., 0., 0.; for (unsigned i = 0; i < solid_acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(solid_acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(liquid_acceleration(i)).epsilon(Tolerance)); } @@ -2208,11 +2208,11 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { true); // Exception check when mass is zero - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nSolid, 0.)); - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nLiquid, 0.)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NSolid, 0.)); + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NLiquid, 0.)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.).epsilon(Tolerance)); REQUIRE(node->compute_acceleration_velocity_twophase_explicit(dt) == false); @@ -2228,90 +2228,90 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { // Check initial momentum for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Check update momentum to 10 REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nSolid, momentum)); + node->update_momentum(true, mpm::NodePhase::NSolid, momentum)); REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nLiquid, momentum)); + node->update_momentum(true, mpm::NodePhase::NLiquid, momentum)); for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Check update momentum to 20 REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nSolid, momentum)); + node->update_momentum(true, mpm::NodePhase::NSolid, momentum)); REQUIRE_NOTHROW( - node->update_momentum(true, mpm::NodePhase::nLiquid, momentum)); + node->update_momentum(true, mpm::NodePhase::NLiquid, momentum)); for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(20.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(20.).epsilon(Tolerance)); } // Check assign momentum to 10 REQUIRE_NOTHROW( - node->update_momentum(false, mpm::NodePhase::nSolid, momentum)); + node->update_momentum(false, mpm::NodePhase::NSolid, momentum)); REQUIRE_NOTHROW( - node->update_momentum(false, mpm::NodePhase::nLiquid, momentum)); + node->update_momentum(false, mpm::NodePhase::NLiquid, momentum)); for (unsigned i = 0; i < momentum.size(); ++i) { - REQUIRE(node->momentum(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NSolid)(i) == Approx(10.).epsilon(Tolerance)); - REQUIRE(node->momentum(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->momentum(mpm::NodePhase::NLiquid)(i) == Approx(10.).epsilon(Tolerance)); } // Check zero velocity for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Check mass double mass = 0.; - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nSolid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NSolid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(0.0).epsilon(Tolerance)); - REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::nLiquid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE_NOTHROW(node->update_mass(false, mpm::NodePhase::NLiquid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(0.0).epsilon(Tolerance)); // Compute and check velocity this should throw zero mass node->compute_velocity(); mass = 100.; // Update mass to 100.5 - REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::nSolid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nSolid) == + REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::NSolid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NSolid) == Approx(100.).epsilon(Tolerance)); - REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::nLiquid, mass)); - REQUIRE(node->mass(mpm::NodePhase::nLiquid) == + REQUIRE_NOTHROW(node->update_mass(true, mpm::NodePhase::NLiquid, mass)); + REQUIRE(node->mass(mpm::NodePhase::NLiquid) == Approx(100.).epsilon(Tolerance)); // Check zero velocity for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } // Compute and check velocity node->compute_velocity(); for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(0.1).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(0.1).epsilon(Tolerance)); } @@ -2320,20 +2320,20 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { for (unsigned i = 0; i < acceleration.size(); ++i) acceleration(i) = 5.; for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(0.).epsilon(Tolerance)); } - REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::nSolid, + REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::NSolid, acceleration)); - REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::nLiquid, + REQUIRE_NOTHROW(node->update_acceleration(true, mpm::NodePhase::NLiquid, acceleration)); for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(5.).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(5.).epsilon(Tolerance)); } @@ -2341,18 +2341,18 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { Eigen::Matrix velocity; velocity << 0.1, 0.1, 0.1; for (unsigned i = 0; i < velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(velocity(i)).epsilon(Tolerance)); } // Check acceleration before constraints acceleration << 5., 5., 5.; for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } SECTION("Check Cartesian velocity constraints") { @@ -2367,17 +2367,17 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { // Check apply constraints velocity << -12.5, 0.1, 0.1; for (unsigned i = 0; i < velocity.size(); ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(velocity(i)).epsilon(Tolerance)); } acceleration << 0., 5., 5.; for (unsigned i = 0; i < acceleration.size(); ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } } @@ -2403,47 +2403,47 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { // Check apply constraints velocity << -14.5068204271, -0.1432759442, 1.4260971922; for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(velocity(i)).epsilon(Tolerance)); } // Check that the velocity is as specified in local coordinate REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nSolid))(0) == + node->velocity(mpm::NodePhase::NSolid))(0) == Approx(-12.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nSolid))(1) == + node->velocity(mpm::NodePhase::NSolid))(1) == Approx(7.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nLiquid))(0) == + node->velocity(mpm::NodePhase::NLiquid))(0) == Approx(-12.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nLiquid))(1) == + node->velocity(mpm::NodePhase::NLiquid))(1) == Approx(7.5).epsilon(Tolerance)); // Check applied constraints on acceleration in the global coordinates acceleration << 0.1998888554, -1.1336260315, 1.9937880031; for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } // Check that the acceleration is 0 in local coordinate REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nSolid))(0) == + node->acceleration(mpm::NodePhase::NSolid))(0) == Approx(0).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nSolid))(1) == + node->acceleration(mpm::NodePhase::NSolid))(1) == Approx(0).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nLiquid))(0) == + node->acceleration(mpm::NodePhase::NLiquid))(0) == Approx(0).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nLiquid))(1) == + node->acceleration(mpm::NodePhase::NLiquid))(1) == Approx(0).epsilon(Tolerance)); } @@ -2471,53 +2471,53 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { // Check apply constraints velocity << 13.351984588153375, -5.717804716697730, 10.572663655835457; for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->velocity(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NSolid)(i) == Approx(velocity(i)).epsilon(Tolerance)); - REQUIRE(node->velocity(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->velocity(mpm::NodePhase::NLiquid)(i) == Approx(velocity(i)).epsilon(Tolerance)); } // Check that the velocity is as specified in local coordinate REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nSolid))(0) == + node->velocity(mpm::NodePhase::NSolid))(0) == Approx(10.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nSolid))(1) == + node->velocity(mpm::NodePhase::NSolid))(1) == Approx(-12.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nSolid))(2) == + node->velocity(mpm::NodePhase::NSolid))(2) == Approx(7.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nLiquid))(0) == + node->velocity(mpm::NodePhase::NLiquid))(0) == Approx(10.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nLiquid))(1) == + node->velocity(mpm::NodePhase::NLiquid))(1) == Approx(-12.5).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->velocity(mpm::NodePhase::nLiquid))(2) == + node->velocity(mpm::NodePhase::NLiquid))(2) == Approx(7.5).epsilon(Tolerance)); // Check apply constraints acceleration << 0, 0, 0; for (unsigned i = 0; i < Dim; ++i) { - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); - REQUIRE(node->acceleration(mpm::NodePhase::nLiquid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NLiquid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } // Check that the acceleration is 0 in local coordinate REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nSolid))(0) == + node->acceleration(mpm::NodePhase::NSolid))(0) == Approx(0).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nSolid))(1) == + node->acceleration(mpm::NodePhase::NSolid))(1) == Approx(0).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nLiquid))(0) == + node->acceleration(mpm::NodePhase::NLiquid))(0) == Approx(0).epsilon(Tolerance)); REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nLiquid))(1) == + node->acceleration(mpm::NodePhase::NLiquid))(1) == Approx(0).epsilon(Tolerance)); } @@ -2533,7 +2533,7 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { // Check apply constraints acceleration << 3.939339828220179, 3.939339828220179, 5.; for (unsigned i = 0; i < acceleration.size(); ++i) - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); } @@ -2556,14 +2556,14 @@ TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") { // Check applied constraints on acceleration in the global coordinates acceleration << 4.602895052828914, 4.492575657560740, 4.751301246937935; for (unsigned i = 0; i < Dim; ++i) - REQUIRE(node->acceleration(mpm::NodePhase::nSolid)(i) == + REQUIRE(node->acceleration(mpm::NodePhase::NSolid)(i) == Approx(acceleration(i)).epsilon(Tolerance)); // Check the acceleration in local coordinate acceleration << 6.878925666702865, 3.365244416454818, 2.302228080558999; for (unsigned i = 0; i < Dim; ++i) REQUIRE((inverse_rotation_matrix * - node->acceleration(mpm::NodePhase::nSolid))(i) == + node->acceleration(mpm::NodePhase::NSolid))(i) == Approx(acceleration(i)).epsilon(Tolerance)); } } diff --git a/tests/particles/particle_twophase_test.cc b/tests/particles/particle_twophase_test.cc index 58d71ee69..024efbe5a 100644 --- a/tests/particles/particle_twophase_test.cc +++ b/tests/particles/particle_twophase_test.cc @@ -913,12 +913,12 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", // Check nodal mass for (unsigned i = 0; i < nodes.size(); ++i) { // Solid phase - REQUIRE(nodes.at(i)->mass(mpm::NodePhase::nSolid) == + REQUIRE(nodes.at(i)->mass(mpm::NodePhase::NSolid) == Approx(nodal_mass.at(i) * (1 - particle->porosity())) .epsilon(Tolerance)); // Liquid phase REQUIRE( - nodes.at(i)->mass(mpm::NodePhase::nLiquid) == + nodes.at(i)->mass(mpm::NodePhase::NLiquid) == Approx(nodal_mass.at(i) * particle->porosity()).epsilon(Tolerance)); } @@ -937,11 +937,11 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", for (unsigned i = 0; i < nodal_momentum.rows(); ++i) for (unsigned j = 0; j < nodal_momentum.cols(); ++j) { // Solid phase - REQUIRE(nodes.at(i)->momentum(mpm::NodePhase::nSolid)(j) == + REQUIRE(nodes.at(i)->momentum(mpm::NodePhase::NSolid)(j) == Approx(nodal_momentum(i, j) * (1 - particle->porosity())) .epsilon(Tolerance)); // Liquid phase - REQUIRE(nodes.at(i)->momentum(mpm::NodePhase::nLiquid)(j) == + REQUIRE(nodes.at(i)->momentum(mpm::NodePhase::NLiquid)(j) == Approx(nodal_momentum(i, j) * particle->porosity()) .epsilon(Tolerance)); } @@ -957,10 +957,10 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", for (unsigned i = 0; i < nodal_velocity.rows(); ++i) for (unsigned j = 0; j < nodal_velocity.cols(); ++j) { // Solid phase - REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::nSolid)(j) == + REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::NSolid)(j) == Approx(nodal_velocity(i, j)).epsilon(Tolerance)); // Liquid phase - REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::nLiquid)(j) == + REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::NLiquid)(j) == Approx(nodal_velocity(i, j)).epsilon(Tolerance)); } @@ -974,11 +974,11 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", for (unsigned i = 0; i < nodes.size(); ++i) { // Solid phase REQUIRE_NOTHROW(nodes.at(i)->update_momentum( - false, mpm::NodePhase::nSolid, + false, mpm::NodePhase::NSolid, nodal_momentum.row(i) * (1 - particle->porosity()))); // Liquid phase REQUIRE_NOTHROW(nodes.at(i)->update_momentum( - false, mpm::NodePhase::nLiquid, + false, mpm::NodePhase::NLiquid, nodal_momentum.row(i) * particle->porosity())); } @@ -995,10 +995,10 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", for (unsigned i = 0; i < nodal_velocity.rows(); ++i) for (unsigned j = 0; j < nodal_velocity.cols(); ++j) { // Solid phase - REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::nSolid)(j) == + REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::NSolid)(j) == Approx(nodal_velocity(i, j)).epsilon(Tolerance)); // Liquid phase - REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::nLiquid)(j) == + REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::NLiquid)(j) == Approx(nodal_velocity(i, j)).epsilon(Tolerance)); } @@ -1069,7 +1069,7 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", // Check nodal body force for (unsigned i = 0; i < body_force.rows(); ++i) for (unsigned j = 0; j < body_force.cols(); ++j) - REQUIRE(nodes[i]->external_force(mpm::NodePhase::nSolid)[j] == + REQUIRE(nodes[i]->external_force(mpm::NodePhase::NSolid)[j] == Approx(body_force(i, j)).epsilon(Tolerance)); // Check traction force @@ -1101,7 +1101,7 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", // Check nodal traction force for (unsigned i = 0; i < traction_force.rows(); ++i) for (unsigned j = 0; j < traction_force.cols(); ++j) - REQUIRE(nodes[i]->external_force(mpm::NodePhase::nSolid)[j] == + REQUIRE(nodes[i]->external_force(mpm::NodePhase::NSolid)[j] == Approx(traction_force(i, j)).epsilon(Tolerance)); // Reset traction particle->assign_traction(direction, @@ -1111,7 +1111,7 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", // Check nodal external force for (unsigned i = 0; i < traction_force.rows(); ++i) for (unsigned j = 0; j < traction_force.cols(); ++j) - REQUIRE(nodes[i]->external_force(mpm::NodePhase::nSolid)[j] == + REQUIRE(nodes[i]->external_force(mpm::NodePhase::NSolid)[j] == Approx(body_force(i, j)).epsilon(Tolerance)); // Internal force @@ -1130,7 +1130,7 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", // Check nodal internal force for (unsigned i = 0; i < internal_force.rows(); ++i) for (unsigned j = 0; j < internal_force.cols(); ++j) - REQUIRE(nodes[i]->internal_force(mpm::NodePhase::nMixture)[j] == + REQUIRE(nodes[i]->internal_force(mpm::NodePhase::NMixture)[j] == Approx(internal_force(i, j)).epsilon(Tolerance)); // Internal force @@ -1170,10 +1170,10 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", for (unsigned i = 0; i < nodal_velocity.rows(); ++i) for (unsigned j = 0; j < nodal_velocity.cols(); ++j) { // Solid phase - REQUIRE(nodes[i]->velocity(mpm::NodePhase::nSolid)[j] == + REQUIRE(nodes[i]->velocity(mpm::NodePhase::NSolid)[j] == Approx(nodal_velocity(i, j)).epsilon(Tolerance)); // Liquid phase - REQUIRE(nodes[i]->velocity(mpm::NodePhase::nLiquid)[j] == + REQUIRE(nodes[i]->velocity(mpm::NodePhase::NLiquid)[j] == Approx(nodal_liquid_velocity(i, j)).epsilon(Tolerance)); } @@ -1194,10 +1194,10 @@ TEST_CASE("TwoPhase Particle is checked for 2D case", for (unsigned i = 0; i < nodal_acceleration.rows(); ++i) for (unsigned j = 0; j < nodal_acceleration.cols(); ++j) { // Solid phase - REQUIRE(nodes[i]->acceleration(mpm::NodePhase::nSolid)[j] == + REQUIRE(nodes[i]->acceleration(mpm::NodePhase::NSolid)[j] == Approx(nodal_acceleration(i, j)).epsilon(Tolerance)); // Liquid phase - REQUIRE(nodes[i]->acceleration(mpm::NodePhase::nLiquid)[j] == + REQUIRE(nodes[i]->acceleration(mpm::NodePhase::NLiquid)[j] == Approx(nodal_liquid_acceleration(i, j)).epsilon(Tolerance)); } // Approx(nodal_velocity(i, j) / dt).epsilon(Tolerance)); @@ -2350,11 +2350,11 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", // Check nodal mass for (unsigned i = 0; i < nodes.size(); ++i) { // Solid phase - REQUIRE(nodes.at(i)->mass(mpm::NodePhase::nSolid) == + REQUIRE(nodes.at(i)->mass(mpm::NodePhase::NSolid) == Approx(nodal_mass.at(i) * (1 - particle->porosity()))); // Liquid phase REQUIRE( - nodes.at(i)->mass(mpm::NodePhase::nLiquid) == + nodes.at(i)->mass(mpm::NodePhase::NLiquid) == Approx(nodal_mass.at(i) * particle->porosity()).epsilon(Tolerance)); } @@ -2378,11 +2378,11 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", for (unsigned i = 0; i < nodal_momentum.rows(); ++i) for (unsigned j = 0; j < nodal_momentum.cols(); ++j) { // Solid phase - REQUIRE(nodes.at(i)->momentum(mpm::NodePhase::nSolid)(j) == + REQUIRE(nodes.at(i)->momentum(mpm::NodePhase::NSolid)(j) == Approx(nodal_momentum(i, j) * (1 - particle->porosity())) .epsilon(Tolerance)); // Liquid phase - REQUIRE(nodes.at(i)->momentum(mpm::NodePhase::nLiquid)(j) == + REQUIRE(nodes.at(i)->momentum(mpm::NodePhase::NLiquid)(j) == Approx(nodal_momentum(i, j) * particle->porosity()) .epsilon(Tolerance)); } @@ -2403,10 +2403,10 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", for (unsigned i = 0; i < nodal_velocity.rows(); ++i) for (unsigned j = 0; j < nodal_velocity.cols(); ++j) { // Solid phase - REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::nSolid)(j) == + REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::NSolid)(j) == Approx(nodal_velocity(i, j)).epsilon(Tolerance)); // Liquid phase - REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::nLiquid)(j) == + REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::NLiquid)(j) == Approx(nodal_velocity(i, j)).epsilon(Tolerance)); } @@ -2424,11 +2424,11 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", for (unsigned i = 0; i < nodes.size(); ++i) { // Solid phase REQUIRE_NOTHROW(nodes.at(i)->update_momentum( - false, mpm::NodePhase::nSolid, + false, mpm::NodePhase::NSolid, nodal_momentum.row(i) * (1 - particle->porosity()))); // Liquid phase REQUIRE_NOTHROW(nodes.at(i)->update_momentum( - false, mpm::NodePhase::nLiquid, + false, mpm::NodePhase::NLiquid, nodal_momentum.row(i) * particle->porosity())); } @@ -2449,10 +2449,10 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", for (unsigned i = 0; i < nodal_velocity.rows(); ++i) for (unsigned j = 0; j < nodal_velocity.cols(); ++j) { // Solid phase - REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::nSolid)(j) == + REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::NSolid)(j) == Approx(nodal_velocity(i, j)).epsilon(Tolerance)); // Liquid phase - REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::nLiquid)(j) == + REQUIRE(nodes.at(i)->velocity(mpm::NodePhase::NLiquid)(j) == Approx(nodal_velocity(i, j)).epsilon(Tolerance)); } @@ -2527,7 +2527,7 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", // Check nodal body force for (unsigned i = 0; i < body_force.rows(); ++i) for (unsigned j = 0; j < body_force.cols(); ++j) - REQUIRE(nodes[i]->external_force(mpm::NodePhase::nSolid)[j] == + REQUIRE(nodes[i]->external_force(mpm::NodePhase::NSolid)[j] == Approx(body_force(i, j)).epsilon(Tolerance)); // Check traction force @@ -2562,7 +2562,7 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", // Check nodal traction force for (unsigned i = 0; i < traction_force.rows(); ++i) for (unsigned j = 0; j < traction_force.cols(); ++j) - REQUIRE(nodes[i]->external_force(mpm::NodePhase::nSolid)[j] == + REQUIRE(nodes[i]->external_force(mpm::NodePhase::NSolid)[j] == Approx(traction_force(i, j)).epsilon(Tolerance)); // Reset traction particle->assign_traction(direction, @@ -2572,7 +2572,7 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", // Check nodal external force for (unsigned i = 0; i < traction_force.rows(); ++i) for (unsigned j = 0; j < traction_force.cols(); ++j) - REQUIRE(nodes[i]->external_force(mpm::NodePhase::nSolid)[j] == + REQUIRE(nodes[i]->external_force(mpm::NodePhase::NSolid)[j] == Approx(body_force(i, j)).epsilon(Tolerance)); // Internal force @@ -2595,12 +2595,12 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", // Check nodal internal force for (unsigned i = 0; i < internal_force.rows(); ++i) for (unsigned j = 0; j < internal_force.cols(); ++j) - REQUIRE(nodes[i]->internal_force(mpm::NodePhase::nMixture)[j] == + REQUIRE(nodes[i]->internal_force(mpm::NodePhase::NMixture)[j] == Approx(internal_force(i, j)).epsilon(Tolerance)); // Calculate nodal acceleration and velocity for (const auto& node : nodes) - node->compute_acceleration_velocity(mpm::NodePhase::nSolid, dt); + node->compute_acceleration_velocity(mpm::NodePhase::NSolid, dt); // Check nodal velocity // clang-format off @@ -2616,7 +2616,7 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", // Check nodal velocity for (unsigned i = 0; i < nodal_velocity.rows(); ++i) for (unsigned j = 0; j < nodal_velocity.cols(); ++j) - REQUIRE(nodes[i]->velocity(mpm::NodePhase::nSolid)[j] == + REQUIRE(nodes[i]->velocity(mpm::NodePhase::NSolid)[j] == Approx(nodal_velocity(i, j)).epsilon(Tolerance)); // Check nodal acceleration @@ -2634,7 +2634,7 @@ TEST_CASE("TwoPhase Particle is checked for 3D case", // Check nodal acceleration for (unsigned i = 0; i < nodal_acceleration.rows(); ++i) for (unsigned j = 0; j < nodal_acceleration.cols(); ++j) - REQUIRE(nodes[i]->acceleration(mpm::NodePhase::nSolid)[j] == + REQUIRE(nodes[i]->acceleration(mpm::NodePhase::NSolid)[j] == Approx(nodal_acceleration(i, j)).epsilon(Tolerance)); // Approx(nodal_velocity(i, j) / dt).epsilon(Tolerance));