diff --git a/include/io/io_mesh.h b/include/io/io_mesh.h index d7e6da99a..b9696e2f2 100644 --- a/include/io/io_mesh.h +++ b/include/io/io_mesh.h @@ -70,6 +70,12 @@ class IOMesh { virtual std::vector> read_particles_scalar_properties(const std::string& scalar_file) = 0; + //! Read particle vector properties + //! \param[in] vector_file file name with particle vector properties + //! \retval Vector of particles vector properties + virtual std::vector>> + read_particles_vector_properties(const std::string& vector_file) = 0; + //! Read pressure constraints file //! \param[in] pressure_constraints_files file name with pressure constraints virtual std::vector> read_pressure_constraints( diff --git a/include/io/io_mesh_ascii.h b/include/io/io_mesh_ascii.h index 9ce82ad82..e115bd58b 100644 --- a/include/io/io_mesh_ascii.h +++ b/include/io/io_mesh_ascii.h @@ -57,6 +57,12 @@ class IOMeshAscii : public IOMesh { std::vector> read_particles_scalar_properties( const std::string& scalar_file) override; + //! Read particle vector properties + //! \param[in] vector_file file name with particle vector properties + //! \retval Vector of particles vector properties + std::vector>> + read_particles_vector_properties(const std::string& vector_file) override; + //! Read pressure constraints file //! \param[in] pressure_constraints_files file name with pressure //! constraints diff --git a/include/io/io_mesh_ascii.tcc b/include/io/io_mesh_ascii.tcc index 303c632ed..f468f5559 100644 --- a/include/io/io_mesh_ascii.tcc +++ b/include/io/io_mesh_ascii.tcc @@ -300,6 +300,55 @@ std::vector> return scalar_properties; } +//! Return particles vector properties +template +std::vector>> + mpm::IOMeshAscii::read_particles_vector_properties( + const std::string& vector_file) { + + // Particles vector properties + std::vector> vector_properties; + + // input file stream + std::fstream file; + file.open(vector_file.c_str(), std::ios::in); + + try { + if (file.is_open() && file.good()) { + // Line + std::string line; + while (std::getline(file, line)) { + boost::algorithm::trim(line); + std::istringstream istream(line); + // ignore comment lines (# or !) or blank lines + if ((line.find('#') == std::string::npos) && + (line.find('!') == std::string::npos) && (line != "")) { + // ID + mpm::Index id; + // Vector + Eigen::Matrix vector; + while (istream.good()) { + // Read stream + istream >> id; + // Read to vector + for (unsigned i = 0; i < Tdim; ++i) istream >> vector[i]; + + vector_properties.emplace_back(std::make_tuple(id, vector)); + } + } + } + } else { + throw std::runtime_error("File not open or not good!"); + } + file.close(); + } catch (std::exception& exception) { + console_->error("Read particle {} #{}: {}\n", __FILE__, __LINE__, + exception.what()); + file.close(); + } + return vector_properties; +} + //! Read pressure constraints file template std::vector> diff --git a/include/mesh/mesh.h b/include/mesh/mesh.h index 59c315d24..f7e6ae370 100644 --- a/include/mesh/mesh.h +++ b/include/mesh/mesh.h @@ -317,10 +317,22 @@ class Mesh { const std::map>& euler_angles); //! Assign particles volumes - //! \param[in] particle_volumes Volume at dir on particle + //! \param[in] particle_volumes Volume on particle bool assign_particles_volumes( const std::vector>& particle_volumes); + //! Assign particles velocities + //! \param[in] particle_velocities Initial particle velocity + bool assign_particles_velocities( + const std::vector>>& + particle_velocities); + + //! Assign particles accelerations + //! \param[in] particle_accelerations Initial particle velocity + bool assign_particles_accelerations( + const std::vector>>& + particle_accelerations); + //! Create particles tractions //! \param[in] mfunction Math function if defined //! \param[in] setid Particle set id diff --git a/include/mesh/mesh.tcc b/include/mesh/mesh.tcc index cdd20d29d..b051fffdc 100644 --- a/include/mesh/mesh.tcc +++ b/include/mesh/mesh.tcc @@ -1313,6 +1313,78 @@ bool mpm::Mesh::assign_particles_volumes( return status; } +//! Assign particle velocities +template +bool mpm::Mesh::assign_particles_velocities( + const std::vector>>& + particle_velocities) { + bool status = true; + + try { + if (!particles_.size()) + throw std::runtime_error( + "No particles have been assigned in mesh, cannot assign particles " + "velocities"); + + if (particles_.size() < particle_velocities.size()) + throw std::runtime_error( + "Number of particles in mesh and initial velocities don't match"); + + // Loop over particle velocities + for (const auto& particle_vel : particle_velocities) { + // Particle id + mpm::Index pid = std::get<0>(particle_vel); + // Velocity vector + VectorDim pvel = std::get<1>(particle_vel); + + if (map_particles_.find(pid) != map_particles_.end()) { + map_particles_[pid]->assign_velocity(pvel); + } + } + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Assign particle accelerations +template +bool mpm::Mesh::assign_particles_accelerations( + const std::vector>>& + particle_accelerations) { + bool status = true; + + try { + if (!particles_.size()) + throw std::runtime_error( + "No particles have been assigned in mesh, cannot assign particles " + "accelerations"); + + if (particles_.size() < particle_accelerations.size()) + throw std::runtime_error( + "Number of particles in mesh and initial accelerations don't match"); + + // Loop over particle accelerations + for (const auto& particle_acc : particle_accelerations) { + // Particle id + mpm::Index pid = std::get<0>(particle_acc); + // Acceleration vector + VectorDim pacc = std::get<1>(particle_acc); + + if (map_particles_.find(pid) != map_particles_.end()) { + map_particles_[pid]->assign_acceleration(pacc); + } + } + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + //! Compute and assign rotation matrix to nodes template bool mpm::Mesh::compute_nodal_rotation_matrices( diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h index 67d26e9fe..6f2f373e4 100644 --- a/include/solvers/mpm_base.h +++ b/include/solvers/mpm_base.h @@ -203,6 +203,20 @@ class MPMBase : public MPM { const Json& mesh_prop, const std::shared_ptr>& particle_io); + //! Particles velocities + //! \param[in] mesh_prop Mesh properties + //! \param[in] particle_io Particle IO handle + void particles_velocities( + const Json& mesh_prop, + const std::shared_ptr>& particle_io); + + //! Particles accelerations + //! \param[in] mesh_prop Mesh properties + //! \param[in] particle_io Particle IO handle + void particles_accelerations( + const Json& mesh_prop, + const std::shared_ptr>& particle_io); + // Particles pore pressures //! \param[in] mesh_prop Mesh properties //! \param[in] particle_io Particle IO handle diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index 8d1fa8ff5..bb35f9505 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -425,6 +425,12 @@ void mpm::MPMBase::initialise_particles() { // Read and assign particles stresses this->particles_stresses(mesh_props, particle_io); + // Read and assign particles initial velocity + this->particles_velocities(mesh_props, particle_io); + + // Read and assign particles initial acceleration + this->particles_accelerations(mesh_props, particle_io); + // Read and assign particles initial pore pressure this->particles_pore_pressures(mesh_props, particle_io); @@ -1444,6 +1450,7 @@ void mpm::MPMBase::particles_volumes( if (mesh_props.find("particles_volumes") != mesh_props.end()) { std::string fparticles_volumes = mesh_props["particles_volumes"].template get(); + if (!io_->file_name(fparticles_volumes).empty()) { bool particles_volumes = mesh_->assign_particles_volumes(particle_io->read_particles_volumes( @@ -1454,12 +1461,117 @@ void mpm::MPMBase::particles_volumes( } } else throw std::runtime_error("Particle volumes JSON data not found"); + } catch (std::exception& exception) { console_->warn("#{}: Particle volumes are undefined; {}", __LINE__, exception.what()); } } +// Read and assign particle velocities +template +void mpm::MPMBase::particles_velocities( + const Json& mesh_props, + const std::shared_ptr>& particle_io) { + try { + // Read particle initial velocities + if (mesh_props.find("particles_velocities") != mesh_props.end()) { + // Get generator type + const std::string type = mesh_props["particles_velocities"]["type"] + .template get(); + if (type == "file") { + std::string fparticle_vel = + mesh_props["particles_velocities"]["location"] + .template get(); + if (!io_->file_name(fparticle_vel).empty()) { + + // Get particle initial velocities + const auto particles_vel = + particle_io->read_particles_vector_properties( + io_->file_name(fparticle_vel)); + + // Assign particle velocities + if (!mesh_->assign_particles_velocities(particles_vel)) + throw std::runtime_error( + "Particles velocities are not properly assigned"); + } + } else if (type == "isotropic") { + Eigen::Matrix in_vel; + in_vel.setZero(); + if (mesh_props["particles_velocities"]["values"].is_array() && + mesh_props["particles_velocities"]["values"].size() == + in_vel.size()) { + for (unsigned i = 0; i < in_vel.size(); ++i) { + in_vel[i] = mesh_props["particles_velocities"]["values"].at(i); + } + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::assign_velocity, + std::placeholders::_1, in_vel)); + } else { + throw std::runtime_error("Initial velocity dimension is invalid"); + } + } + } else + throw std::runtime_error("Particle velocities JSON data not found"); + + } catch (std::exception& exception) { + console_->warn("#{}: Particle velocities are undefined {} ", __LINE__, + exception.what()); + } +} + +// Read and assign particle accelerations +template +void mpm::MPMBase::particles_accelerations( + const Json& mesh_props, + const std::shared_ptr>& particle_io) { + try { + // Read particle initial accelerations + if (mesh_props.find("particles_accelerations") != mesh_props.end()) { + // Get generator type + const std::string type = mesh_props["particles_accelerations"]["type"] + .template get(); + if (type == "file") { + std::string fparticle_acc = + mesh_props["particles_accelerations"]["location"] + .template get(); + if (!io_->file_name(fparticle_acc).empty()) { + + // Get particle initial accelerations + const auto particles_acc = + particle_io->read_particles_vector_properties( + io_->file_name(fparticle_acc)); + + // Assign particle accelerations + if (!mesh_->assign_particles_accelerations(particles_acc)) + throw std::runtime_error( + "Particles accelerations are not properly assigned"); + } + } else if (type == "isotropic") { + Eigen::Matrix in_acc; + in_acc.setZero(); + if (mesh_props["particles_accelerations"]["values"].is_array() && + mesh_props["particles_accelerations"]["values"].size() == + in_acc.size()) { + for (unsigned i = 0; i < in_acc.size(); ++i) { + in_acc[i] = mesh_props["particles_accelerations"]["values"].at(i); + } + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::assign_acceleration, + std::placeholders::_1, in_acc)); + } else { + throw std::runtime_error("Initial acceleration dimension is invalid"); + } + } + } else + throw std::runtime_error("Particle accelerations JSON data not found"); + + } catch (std::exception& exception) { + console_->warn("#{}: Particle accelerations are undefined {} ", __LINE__, + exception.what()); + } +} + // Particle velocity constraints template void mpm::MPMBase::particle_velocity_constraints() { diff --git a/tests/io/io_mesh_ascii_test.cc b/tests/io/io_mesh_ascii_test.cc index 58356a53c..23dadac1d 100644 --- a/tests/io/io_mesh_ascii_test.cc +++ b/tests/io/io_mesh_ascii_test.cc @@ -863,6 +863,59 @@ TEST_CASE("IOMeshAscii is checked for 2D", "[IOMesh][IOMeshAscii][2D]") { } } + SECTION("Check particles vector properties file") { + // Particle vector properties + std::map> particles_vectors; + particles_vectors.emplace( + std::make_pair(0, (Eigen::Matrix(10.5, 20.5)))); + particles_vectors.emplace( + std::make_pair(1, (Eigen::Matrix(30.5, -40.5)))); + particles_vectors.emplace( + std::make_pair(2, (Eigen::Matrix(-50.5, -60.5)))); + particles_vectors.emplace( + std::make_pair(3, (Eigen::Matrix(-70.5, 80.5)))); + + // Dump particle vector properties as an input file to be read + std::ofstream file; + file.open("particles-vectors-2d.txt"); + // Write particle vector properties + for (const auto& vectors : particles_vectors) { + file << vectors.first << "\t"; + for (unsigned i = 0; i < dim; ++i) file << (vectors.second)(i) << "\t"; + file << "\n"; + } + + file.close(); + + // Check read particles vector properties file + SECTION("Check read particles vector properties file") { + // Create a io_mesh object + auto io_mesh = std::make_unique>(); + + // Try to read particles vector properties from a non-existant file + auto read_particles_vectors = io_mesh->read_particles_vector_properties( + "particles-vectors-missing.txt"); + + // Check number of particles vector properties + REQUIRE(read_particles_vectors.size() == 0); + + // Check particles vector properties + read_particles_vectors = + io_mesh->read_particles_vector_properties("particles-vectors-2d.txt"); + + // Check number of particles + REQUIRE(read_particles_vectors.size() == particles_vectors.size()); + + // Check particles vector properties + for (unsigned i = 0; i < particles_vectors.size(); ++i) { + for (unsigned j = 0; j < particles_vectors.at(i).size(); ++j) { + REQUIRE(std::get<1>(read_particles_vectors.at(i))[j] == + Approx(particles_vectors.at(i)[j]).epsilon(Tolerance)); + } + } + } + } + SECTION("Check math function file") { // Vector of math function entries std::array, 2> entries; @@ -1812,6 +1865,59 @@ TEST_CASE("IOMeshAscii is checked for 3D", "[IOMesh][IOMeshAscii][3D]") { } } + SECTION("Check particles vector properties file") { + // Particle vector properties + std::map> particles_vectors; + particles_vectors.emplace( + std::make_pair(0, (Eigen::Matrix(10.5, 20.5, 30.5)))); + particles_vectors.emplace( + std::make_pair(1, (Eigen::Matrix(40.5, -50.5, -60.5)))); + particles_vectors.emplace(std::make_pair( + 2, (Eigen::Matrix(-70.5, -80.5, -90.5)))); + particles_vectors.emplace(std::make_pair( + 3, (Eigen::Matrix(-100.5, 110.5, 120.5)))); + + // Dump particle vector properties as an input file to be read + std::ofstream file; + file.open("particles-vectors-3d.txt"); + // Write particle vector properties + for (const auto& vectors : particles_vectors) { + file << vectors.first << "\t"; + for (unsigned i = 0; i < vectors.second.size(); ++i) + file << (vectors.second)(i) << "\t"; + file << "\n"; + } + + file.close(); + + // Check read particles vector properties file + SECTION("Check read particles vector properties file") { + // Create a io_mesh object + auto io_mesh = std::make_unique>(); + + // Try to read particles vector properties from a non-existant file + auto read_particles_vectors = io_mesh->read_particles_vector_properties( + "particles-vector-missing.txt"); + // Check number of particles vector properties + REQUIRE(read_particles_vectors.size() == 0); + + // Check particles vector properties + read_particles_vectors = + io_mesh->read_particles_vector_properties("particles-vectors-3d.txt"); + + // Check number of particles + REQUIRE(read_particles_vectors.size() == particles_vectors.size()); + + // Check particles vector properties + for (unsigned i = 0; i < particles_vectors.size(); ++i) { + for (unsigned j = 0; j < particles_vectors.at(i).size(); ++j) { + REQUIRE(std::get<1>(read_particles_vectors.at(i))[j] == + Approx(particles_vectors.at(i)[j]).epsilon(Tolerance)); + } + } + } + } + SECTION("Check math function file") { // Vector of math function entries std::array, 2> entries; diff --git a/tests/io/write_mesh_particles_unitcell.cc b/tests/io/write_mesh_particles_unitcell.cc index c6ffe32ae..50817501c 100644 --- a/tests/io/write_mesh_particles_unitcell.cc +++ b/tests/io/write_mesh_particles_unitcell.cc @@ -242,6 +242,8 @@ bool write_json_unitcell_implicit(unsigned dim, const std::string& analysis, std::vector xvalues{{0.0, 0.5, 1.0}}; std::vector fxvalues{{0.0, 1.0, 1.0}}; std::vector initial_stress{{-1.0, -1.0, -1.0, 0.0}}; + std::vector initial_vel{{1.0, 0.}}; + std::vector initial_acc{{0.5, 0.}}; // 3D if (dim == 3) { @@ -255,6 +257,8 @@ bool write_json_unitcell_implicit(unsigned dim, const std::string& analysis, gravity.clear(); gravity = {0., 0., -9.81}; initial_stress = {{-1.0, -1.0, -1.0, 0.0, 0.0, 0.0}}; + initial_vel = {{1.0, 0., 0.}}; + initial_acc = {{0.5, 0., 0.}}; } Json json_file = { @@ -271,7 +275,12 @@ bool write_json_unitcell_implicit(unsigned dim, const std::string& analysis, {"nonlocal_mesh_properties", {{"type", "LME"}, {"gamma", 3}, {"anisotropy", true}}}, {"particle_stresses", - {{"type", "isotropic"}, {"value", initial_stress}}}}}, + {{"type", "isotropic"}, {"value", initial_stress}}}, + {"particle_velocities", + {{"type", "isotropic"}, {"value", initial_vel}}}, + {"particle_accelerations", + {{"type", "isotropic"}, {"value", initial_acc}}} + }}, {"particles", {{{"generator", {{"type", "file"}, diff --git a/tests/mesh/mesh_test_2d.cc b/tests/mesh/mesh_test_2d.cc index 9a7512462..29850886d 100644 --- a/tests/mesh/mesh_test_2d.cc +++ b/tests/mesh/mesh_test_2d.cc @@ -1119,6 +1119,79 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { particles_stresses) == false); } + // Test assign initial particles velocities + SECTION("Check assign initial particles velocities") { + // Vector of particle velocities + std::vector>> + particles_velocities; + + // Check number of particles in mesh + REQUIRE(mesh->nparticles() == 8); + + // Velocities + particles_velocities.emplace_back(std::make_tuple( + 0, Eigen::Matrix::Constant(0.0))); + particles_velocities.emplace_back(std::make_tuple( + 1, Eigen::Matrix::Constant(0.1))); + particles_velocities.emplace_back(std::make_tuple( + 2, Eigen::Matrix::Constant(-0.2))); + particles_velocities.emplace_back(std::make_tuple( + 3, Eigen::Matrix::Constant(0.3))); + particles_velocities.emplace_back(std::make_tuple( + 4, Eigen::Matrix::Constant(-0.4))); + particles_velocities.emplace_back(std::make_tuple( + 5, Eigen::Matrix::Constant(0.5))); + particles_velocities.emplace_back(std::make_tuple( + 6, Eigen::Matrix::Constant(-0.6))); + particles_velocities.emplace_back(std::make_tuple( + 7, Eigen::Matrix::Constant(-0.7))); + + REQUIRE(mesh->assign_particles_velocities(particles_velocities) == + true); + + // When velocities fail + particles_velocities.emplace_back(std::make_tuple( + 8, Eigen::Matrix::Constant(0.8))); + REQUIRE(mesh->assign_particles_velocities(particles_velocities) == + false); + } + + // Test assign initial particles accelerations + SECTION("Check assign initial particles accelerations") { + // Vector of particle accelerations + std::vector>> + particles_accelerations; + + // Check number of particles in mesh + REQUIRE(mesh->nparticles() == 8); + // Accelerations + particles_accelerations.emplace_back(std::make_tuple( + 0, Eigen::Matrix::Constant(0.0))); + particles_accelerations.emplace_back(std::make_tuple( + 1, Eigen::Matrix::Constant(0.1))); + particles_accelerations.emplace_back(std::make_tuple( + 2, Eigen::Matrix::Constant(-0.2))); + particles_accelerations.emplace_back(std::make_tuple( + 3, Eigen::Matrix::Constant(0.3))); + particles_accelerations.emplace_back(std::make_tuple( + 4, Eigen::Matrix::Constant(-0.4))); + particles_accelerations.emplace_back(std::make_tuple( + 5, Eigen::Matrix::Constant(0.5))); + particles_accelerations.emplace_back(std::make_tuple( + 6, Eigen::Matrix::Constant(-0.6))); + particles_accelerations.emplace_back(std::make_tuple( + 7, Eigen::Matrix::Constant(-0.7))); + + REQUIRE(mesh->assign_particles_accelerations( + particles_accelerations) == true); + + // When accelerations fail + particles_accelerations.emplace_back(std::make_tuple( + 8, Eigen::Matrix::Constant(0.8))); + REQUIRE(mesh->assign_particles_accelerations( + particles_accelerations) == false); + } + // Test assign particles velocity constraints SECTION("Check assign particles velocity constraints") { tsl::robin_map> particle_sets;