diff --git a/README.md b/README.md index 89f77f3..7645b3f 100644 --- a/README.md +++ b/README.md @@ -58,9 +58,9 @@ With _Matter_, you can simulate granular flow on various simple and complex terr * Supports no-slip as well as frictional **boundary conditions** which can be supplied a **Coulomb friction parameter** - * No-slip (`NOSLIP`) - * Frictional slip free to separate (`SLIPFREE`) - * Frictional slip constrained to boundary (`SLIPSTICK`) + * No-slip (`BC::NoSlip`) + * Frictional slip free to separate (`BC::SlipFree`) + * Frictional slip constrained to boundary (`BC::SlipStick`) * **Material-induced boundary friction (MIBF)**: A potentially variable internal friction parameter of a particular plastic model is used as the Coulomb friction for terrain-material interaction @@ -109,9 +109,13 @@ Download [VSCode](https://code.visualstudio.com/) and [Docker](https://www.docke 5. Compile (NB: the number of threads for the _simulation_ is specified in `mpm.cpp`) `make -j ` -6. Run the executable: +6. Run the sim you set up in `mpm.cpp`: `./src/mpm` +7. [Optional] Tests: + Run all tests with `make test` or `ctest`. Run single test with `ctest -R ` + + ### Example of setup file Here is a minimal setup file `mpm.cpp` for a simple granular collapse. @@ -128,27 +132,32 @@ int main(){ sim.initialize(/*save to file*/ true, /*path*/ "output/", /*name*/ "collapse"); - sim.end_frame = 20; // last frame to simulate - sim.fps = 10; // frames per second - sim.n_threads = 8; // number of threads in parallel + sim.end_frame = 20; // last frame to simulate + sim.fps = 10; // frames per second + sim.n_threads = 8; // number of threads in parallel sim.Lx = 1; sim.Ly = 1; - sim.Lz = 1; // ONLY IF 3D, OTHERWISE REMOVE LINE - SampleParticles(sim, /*sampling radius*/ 0.01); + // sim.Lz = 1; // only 3D, otherwise remove + + sampleParticles(sim, /*sampling radius*/ 0.01); - ObjectPlate ground = ObjectPlate(/*position*/ 0, /*type*/ bottom, /*boundary cond.*/ SLIPFREE, /*friction*/ 0.5); - sim.plates.push_back(ground); + sim.plates.push_back(std::make_unique( + /*position*/ 0, + /*type*/ PlateType::bottom, + /*BC*/ BC::NoSlip, + /*friction*/ 0.5 + )); - sim.rho = 1000; // Density (kg/m3) - sim.gravity[1] = -9.81; // Gravity + sim.rho = 1000; // density (kg/m3) + sim.gravity[1] = -9.81; // gravity sim.E = 1e6; // Young's modulus (Pa) sim.nu = 0.3; // Poisson's ratio (-) - sim.plastic_model = DP; // Drucker-Prager yield - sim.M = 0.5; // Internal friction - sim.q_cohesion = 0; // Cohesion + sim.plastic_model = PlasticModel::DP; // Drucker-Prager yield + sim.M = 0.5; // internal friction + sim.q_cohesion = 0; // cohesion sim.simulate(); return 0; @@ -164,7 +173,7 @@ Rigid objects and terrains (boundaries) are either * or imported as `.vdb` level sets files using [OpenVDB](https://www.openvdb.org/) ##### Analytical objects -Analytical objects can be specified as a derived class from the general `ObjectGeneral` class. An example of this is `ObjectBump` which provides the terrain of a smooth bump used in the flow experiments in [Viroulet et al. (2017)](https://doi.org/10.1017/jfm.2017.41). For the very common case of an axis-aligned plate, an `ObjectPlate` class has been made separate from `ObjectGeneral` class for efficiency and convenience. In `ObjectPlate`, you can also assign a speed to the plate, as well as controls on the time-evolution of the speed. Any plate must either a `top`, `bottom`, `front`, `back`, `left` or `right`. +Analytical objects can be specified as a derived class from the general `ObjectGeneral` class. An example of this is `ObjectBump` which provides the terrain of a smooth bump used in the flow experiments in [Viroulet et al. (2017)](https://doi.org/10.1017/jfm.2017.41). For the very common case of an axis-aligned plate, an `ObjectPlate` class has been made separate from `ObjectGeneral` class for efficiency and convenience. In `ObjectPlate`, you can also assign a speed to the plate, as well as controls on the time-evolution of the speed. Any plate must either a `PlateType::top`, `PlateType::bottom`, `PlateType::front`, `PlateType::back`, `PlateType::left` or `PlateType::right`. ##### OpenVDB objects A terrain/object from a `.vdb` is stored in an instance of the `ObjectVdb` class which is derived from `ObjectGeneral`. @@ -199,69 +208,69 @@ This is a non-exhaustive list of parameters and options (of the `Simulation` cla | `reduce_verbose` | false | Reduce writing to screen | `use_musl` | false | Use MUSL instead of USL | `use_mibf` | false | Use Material-Induced Boundary Friction (MIBF), only relevant for certain plasticity models -| `use_mises_q` | false | If `true` define the "equivalent shear stress" q as the von Mises equivalent stress q = sqrt(3/2 s:s), otherwise q = sqrt(1/2 s:s). | `pbc` | false | Use periodic boundary conditions in $x$-direction bounded by `Lx` -| `Lx`, `Ly`, `Lz` | 1.0 | The material sample space used in `SampleParticles(...)` +| `Lx`, `Ly`, `Lz` | 1.0 | The material sample space used in `sampleParticles(...)`. Not needed when sampling from VDB. | `grid_reference_point` | - | Optionally provide a point to be considered in the initial adaptive grid creation, otherwise it only considers the particle domain -| `elastic_model` | Hencky | Elastic model. Note that Hencky's model must be used when combined with a plastic model. -| `plastic_model` | NoPlasticity | Plastic model. Parameters are set according to the model used, see below. +| `elastic_model` | ElasticModel::Hencky | Elastic model. Note that Hencky's model must be used when combined with a plastic model. +| `plastic_model` | PlasticModel::NoPlasticity | Plastic model. Parameters are set according to the model used, see below. | `E` | 1e6 | The 3D Young's modulus (Pa) | `nu` | 0.3 | The 3D Poisson's ratio (-) +| `q_prefac` | sqrt(1/2) | Prefactor used in definition of q, default is q = sqrt(1/2 s:s). Here is a list of the various plastic models and their parameters: | Model | Name | Parameters | Default value | | ---- | ---- | ---- | --- | -| Von Mises | `VM` | `q_max` | 100.0 | -| Drucker-Prager | `DP` | `M` | 1.0 | +| Von Mises | `PlasticModel::VM` | `q_max` | 100.0 | +| Drucker-Prager | `PlasticModel::DP` | `M` | 1.0 | | | | `q_cohesion` | 0.0 | -| Drucker-Prager with strain-softening | `DPSoft` | `M` | 1.0 | +| Drucker-Prager strain-softening | `PlasticModel::DPSoft` | `M` | 1.0 | | | | `q_cohesion` | 0.0 | | | | `xi` | 0.0 | | | | `use_pradhana` | true | -| Modified Cam-Clay | `MCC` | `beta` | 0.0 | +| Modified Cam-Clay | `PlasticModel::MCC` | `beta` | 0.0 | | | | `p0` | 1000.0 | | | | `xi` | 0.0 | | | | `M` | 1.0 | -| Perzyna-Von Mises | `VMVisc` | `q_max` | 100.0 | +| Perzyna-Von Mises | `PlasticModel::VMVisc` | `q_max` | 100.0 | | | | `q_min` | 100.0 | | | | `p_min` | -1.0e20 | | | | `xi` | 0.0 | | | | `perzyna_exp` | 1.0 | | | | `perzyna_visc` | 0.0 | -| Perzyna-Drucker-Prager | `DPVisc` | `M` | 1.0 | +| Perzyna-Drucker-Prager | `PlasticModel::DPVisc` | `M` | 1.0 | | | | `q_cohesion` | 0.0 | | | | `use_pradhana` | true | | | | `perzyna_exp` | 1.0 | | | | `perzyna_visc` | 0.0 | -| Perzyna-Modified Cam-Clay | `MCCVisc` | `beta` | 0.0 | +| Perzyna-Modified Cam-Clay | `PlasticModel::MCCVisc` | `beta` | 0.0 | | | | `p0` | 1000.0 | | | | `xi` | 0.0 | | | | `M` | 1.0 | | | | `perzyna_exp` | 1.0 | | | | `perzyna_visc` | 0.0 | -| $\mu(I)$-rheology | `DPMui` | `q_cohesion` | 0.0 | +| $\mu(I)$-rheology | `PlasticModel::DPMui` | `q_cohesion` | 0.0 | | | | `use_pradhana` | true | -| | | `rho_s` | 1.0 | +| | | `rho_s` | 2500 | | | | `grain_diameter`| 0.001 | | | | `I_ref` | 0.279 | | | | `mu_1` | 0.382 | | | | `mu_2` | 0.644 | -| Critical state $\mu(I)$-rheology | `MCCMui` | `beta` | 0.0 | +| Critical state $\mu(I)$-rheology | `PlasticModel::MCCMui` | `beta` | 0.0 | | | | `p0` | 1000.0 | | | | `xi` | 0.0 | -| | | `rho_s` | 1000.0 | +| | | `rho_s` | 2500 | | | | `grain_diameter`| 0.001 | | | | `I_ref` | 0.279 | | | | `mu_1` | 0.382 | | | | `mu_2` | 0.644 | In the MCC-based models, one must also choose a corresponding hardening law: -* exponential explicit law `ExpoExpl` -* hyperbolic sine explicit law `SinhExpl` -* exponential implicit law `ExpoImpl` -* hyperbolic sine implicit law `SinhImpl` +* exponential explicit law `HardeningLaw::ExpoExpl` +* hyperbolic sine explicit law `HardeningLaw::SinhExpl` +* exponential implicit law `HardeningLaw::ExpoImpl` +* hyperbolic sine implicit law `HardeningLaw::SinhImpl` ## Limitations @@ -280,7 +289,9 @@ In the MCC-based models, one must also choose a corresponding hardening law: ## License and attribution _Matter_ is an open-source software licensed under _GNU General Public License v3.0_ (see LICENSE file). -If you are interested in using Matter in commercial products or services, please do not hesitate to contact [Lars Blatny](https://larsblatny.github.io/) (lars.blatny [at] slf.ch). +If you are interested in using Matter in commercial products or services, please do not hesitate to contact [Lars Blatny](https://larsblatny.github.io/) (lars.blatny [at] slf.ch). + +Please attribute this software by citing this article: [doi:10.5194/egusphere-2025-1157](https://doi.org/10.5194/egusphere-2025-1157) -If you use the $\mu(I)$-rheology models, please cite this article: [doi:10.1017/jfm.2024.643](https://doi.org/10.1017/jfm.2024.643) +If you use the $\mu(I)$-rheology models, please cite this article: [doi:10.1017/jfm.2024.643](https://doi.org/10.1017/jfm.2024.643) diff --git a/examples/collapse.cpp b/examples/collapse.cpp index d2004ef..8c753d2 100644 --- a/examples/collapse.cpp +++ b/examples/collapse.cpp @@ -29,7 +29,7 @@ int main(){ sim.flip_ratio = -0.95; // (A)PIC-(A)FLIP ratio in [-1,1]. // INITILIZE ELASTICITY - sim.elastic_model = Hencky; + sim.elastic_model = ElasticModel::Hencky; sim.E = 1e6; // Young's modulus (Pa) sim.nu = 0.3; // Poisson's ratio (-) sim.rho = 1000; // Density (kg/m3) @@ -61,21 +61,21 @@ int main(){ // sim.particles.v = ... ////// OBJECTS AND TERRAINS - ObjectPlate ground = ObjectPlate(0, bottom, NOSLIP); sim.plates.push_back(ground); + sim.plates.push_back(std::make_unique(0, PlateType::bottom, BC::NoSlip)); /////// Here are some examples how to use the objects derived from ObjectGeneral: // T friction = 0.2; - // ObjectBump bump = ObjectBump(SLIPFREE, friction); sim.objects.push_back(&bump); - // ObjectGate gate = ObjectGate(SLIPFREE, friction); sim.objects.push_back(&gate); + // sim.objects.push_back(std::make_unique(BC::SlipFree, friction)); + // sim.objects.push_back(std::make_unique(BC::SlipFree, friction)); - /////// Here is an example how to use ObjectVdb (uncomment include files and openvdb::initialize() function above): - // ObjectVdb terrain = ObjectVdb("../levelsets/vdb_file_name.vdb", NOSLIP, friction); sim.objects.push_back(&terrain); + /////// Here is an example how to use ObjectVdb (uncomment includes and openvdb::initialize() above): + // sim.objects.push_back(std::make_unique("../levelsets/vdb_file_name.vdb", BC::NoSlip, friction)); ////// PLASTICITY - sim.plastic_model = DPVisc; // Perzyna model with Drucker_Prager yield surface + sim.plastic_model = PlasticModel::DPVisc; // Perzyna model with Drucker_Prager yield surface sim.use_pradhana = true; // Supress unwanted volume expansion in Drucker-Prager models - sim.use_mises_q = false; // [default: false] if true, q is defined as q = sqrt(3/2 * s:s), otherwise q = sqrt(1/2 * s:s) + sim.q_prefac = 1.0 / std::sqrt(2.0); // [default: sqrt(1/2)] Prefactor in def. of q, here q = sqrt(1/2 * s:s) sim.M = std::tan(30*M_PI/180.0); // Internal friction sim.q_cohesion = 0; // Yield surface's intercection of q-axis (in Pa), 0 is the cohesionless case diff --git a/examples/cube_rotating.cpp b/examples/cube_rotating.cpp index 7723928..cb72c5c 100755 --- a/examples/cube_rotating.cpp +++ b/examples/cube_rotating.cpp @@ -1,8 +1,8 @@ // Copyright (C) 2024 Lars Blatny. Released under GPL-3.0 license. -#include "simulation.hpp" #include "tools.hpp" -#include "sampling_particles.hpp" +#include "simulation/simulation.hpp" +#include "sampling/sampling_particles.hpp" int main(){ @@ -43,8 +43,8 @@ int main(){ } // Elasticity - sim.elastic_model = Hencky; - sim.plastic_model = NoPlasticity; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::NoPlasticity; sim.E = 1e6; // Young's modulus (Pa) sim.nu = 0.3; // Poisson's ratio (-) sim.rho = 1550; // Density (kg/m3) diff --git a/src/data_structures.hpp b/src/data_structures.hpp index 5b873c5..16b045e 100644 --- a/src/data_structures.hpp +++ b/src/data_structures.hpp @@ -8,22 +8,20 @@ class Particles{ public: Particles(unsigned int Np = 1){ - x.resize(Np); std::fill( x.begin(), x.end(), TV::Zero() ); - v.resize(Np); std::fill( v.begin(), v.end(), TV::Zero() ); - pic.resize(Np); std::fill( pic.begin(), pic.end(), TV::Zero() ); + x.resize(Np); std::fill( x.begin(), x.end(), TV::Zero() ); + v.resize(Np); std::fill( v.begin(), v.end(), TV::Zero() ); + pic.resize(Np); std::fill( pic.begin(), pic.end(), TV::Zero() ); flip.resize(Np); std::fill( flip.begin(), flip.end(), TV::Zero() ); - eps_pl_dev.resize(Np); std::fill( eps_pl_dev.begin(), eps_pl_dev.end(), 0.0 ); - eps_pl_vol.resize(Np); std::fill( eps_pl_vol.begin(), eps_pl_vol.end(), 0.0 ); - eps_pl_vol_pradhana.resize(Np);std::fill( eps_pl_vol_pradhana.begin(),eps_pl_vol_pradhana.end(),0.0 ); - - delta_gamma.resize(Np); std::fill( delta_gamma.begin(), delta_gamma.end(), 0.0 ); + eps_pl_dev.resize(Np); std::fill( eps_pl_dev.begin(), eps_pl_dev.end(), 0.0 ); + eps_pl_vol.resize(Np); std::fill( eps_pl_vol.begin(), eps_pl_vol.end(), 0.0 ); + eps_pl_vol_pradhana.resize(Np);std::fill( eps_pl_vol_pradhana.begin(),eps_pl_vol_pradhana.end(), 0.0 ); + delta_gamma.resize(Np); std::fill( delta_gamma.begin(), delta_gamma.end(), 0.0 ); viscosity.resize(Np); std::fill( viscosity.begin(), viscosity.end(), 0.0 ); - muI.resize(Np); std::fill( muI.begin(), muI.end(), 0.0 ); + muI.resize(Np); std::fill( muI.begin(), muI.end(), 0.0 ); - tau.resize(Np); std::fill( tau.begin(), tau.end(), TM::Zero() ); - F.resize(Np); std::fill( F.begin(), F.end(), TM::Identity() ); - Bmat.resize(Np); std::fill( Bmat.begin(), Bmat.end(), TM::Zero() ); + F.resize(Np); std::fill( F.begin(), F.end(), TM::Identity() ); + Bmat.resize(Np); std::fill( Bmat.begin(), Bmat.end(), TM::Zero() ); } std::vector x; @@ -34,12 +32,10 @@ class Particles{ std::vector eps_pl_dev; std::vector eps_pl_vol; std::vector eps_pl_vol_pradhana; - std::vector delta_gamma; std::vector viscosity; std::vector muI; - std::vector tau; std::vector F; std::vector Bmat; diff --git a/src/mpm.cpp b/src/mpm.cpp index d2004ef..8c753d2 100644 --- a/src/mpm.cpp +++ b/src/mpm.cpp @@ -29,7 +29,7 @@ int main(){ sim.flip_ratio = -0.95; // (A)PIC-(A)FLIP ratio in [-1,1]. // INITILIZE ELASTICITY - sim.elastic_model = Hencky; + sim.elastic_model = ElasticModel::Hencky; sim.E = 1e6; // Young's modulus (Pa) sim.nu = 0.3; // Poisson's ratio (-) sim.rho = 1000; // Density (kg/m3) @@ -61,21 +61,21 @@ int main(){ // sim.particles.v = ... ////// OBJECTS AND TERRAINS - ObjectPlate ground = ObjectPlate(0, bottom, NOSLIP); sim.plates.push_back(ground); + sim.plates.push_back(std::make_unique(0, PlateType::bottom, BC::NoSlip)); /////// Here are some examples how to use the objects derived from ObjectGeneral: // T friction = 0.2; - // ObjectBump bump = ObjectBump(SLIPFREE, friction); sim.objects.push_back(&bump); - // ObjectGate gate = ObjectGate(SLIPFREE, friction); sim.objects.push_back(&gate); + // sim.objects.push_back(std::make_unique(BC::SlipFree, friction)); + // sim.objects.push_back(std::make_unique(BC::SlipFree, friction)); - /////// Here is an example how to use ObjectVdb (uncomment include files and openvdb::initialize() function above): - // ObjectVdb terrain = ObjectVdb("../levelsets/vdb_file_name.vdb", NOSLIP, friction); sim.objects.push_back(&terrain); + /////// Here is an example how to use ObjectVdb (uncomment includes and openvdb::initialize() above): + // sim.objects.push_back(std::make_unique("../levelsets/vdb_file_name.vdb", BC::NoSlip, friction)); ////// PLASTICITY - sim.plastic_model = DPVisc; // Perzyna model with Drucker_Prager yield surface + sim.plastic_model = PlasticModel::DPVisc; // Perzyna model with Drucker_Prager yield surface sim.use_pradhana = true; // Supress unwanted volume expansion in Drucker-Prager models - sim.use_mises_q = false; // [default: false] if true, q is defined as q = sqrt(3/2 * s:s), otherwise q = sqrt(1/2 * s:s) + sim.q_prefac = 1.0 / std::sqrt(2.0); // [default: sqrt(1/2)] Prefactor in def. of q, here q = sqrt(1/2 * s:s) sim.M = std::tan(30*M_PI/180.0); // Internal friction sim.q_cohesion = 0; // Yield surface's intercection of q-axis (in Pa), 0 is the cohesionless case diff --git a/src/objects/object_bump.hpp b/src/objects/object_bump.hpp index 0141bd7..7c878db 100644 --- a/src/objects/object_bump.hpp +++ b/src/objects/object_bump.hpp @@ -10,9 +10,9 @@ class ObjectBump : public ObjectGeneral { ~ObjectBump(){} - ObjectBump(BoundaryCondition bc_in, T friction_in, std::string name_in = "") : ObjectGeneral(bc_in, friction_in, name_in) {} + ObjectBump(BC bc_in, T friction_in, std::string name_in = "") : ObjectGeneral(bc_in, friction_in, name_in) {} - bool inside(const TV& X_in) override { + bool inside(const TV& X_in) const override { T x = X_in(0); T y = X_in(1); @@ -26,10 +26,10 @@ class ObjectBump : public ObjectGeneral { } - TV normal(const TV& X_in) override { + TV normal(const TV& X_in) const override { T x = X_in(0); - TV n; + TV n = TV::Zero(); T arg = 25*(x-0.43); T b_der = -25 * 0.0475 * std::tanh(arg) / std::cosh(arg); diff --git a/src/objects/object_curve.hpp b/src/objects/object_curve.hpp index 07cc861..18b6b67 100644 --- a/src/objects/object_curve.hpp +++ b/src/objects/object_curve.hpp @@ -10,9 +10,9 @@ class ObjectCurve : public ObjectGeneral { ~ObjectCurve(){} - ObjectCurve(BoundaryCondition bc_in = NOSLIP, T friction_in = 0.0, std::string name_in = "") : ObjectGeneral(bc_in, friction_in, name_in) {} + ObjectCurve(BC bc_in = BC::NoSlip, T friction_in = 0.0, std::string name_in = "") : ObjectGeneral(bc_in, friction_in, name_in) {} - bool inside(const TV& X_in) override { + bool inside(const TV& X_in) const override { T x = X_in(0); T y = X_in(1); @@ -24,10 +24,10 @@ class ObjectCurve : public ObjectGeneral { } - TV normal(const TV& X_in) override { + TV normal(const TV& X_in) const override { T x = X_in(0); - TV n; + TV n = TV::Zero(); n(0) = -2*x; n(1) = 1; diff --git a/src/objects/object_gate.hpp b/src/objects/object_gate.hpp index e74ad1f..9567290 100644 --- a/src/objects/object_gate.hpp +++ b/src/objects/object_gate.hpp @@ -11,9 +11,9 @@ class ObjectGate : public ObjectGeneral { ~ObjectGate(){} - ObjectGate(BoundaryCondition bc_in, T friction_in, std::string name_in = "", T height_in = 0.016) : ObjectGeneral(bc_in, friction_in, name_in), height(height_in) {} + ObjectGate(BC bc_in, T friction_in, std::string name_in = "", T height_in = 0.016) : ObjectGeneral(bc_in, friction_in, name_in), height(height_in) {} - bool inside(const TV& X_in) override { + bool inside(const TV& X_in) const override { T x = X_in(0); T y = X_in(1); @@ -27,10 +27,10 @@ class ObjectGate : public ObjectGeneral { } - TV normal(const TV& X_in) override { + TV normal(const TV& X_in) const override { T x = X_in(0); - TV n; + TV n = TV::Zero(); T b_der = 2.0 * 100 * x; n(0) = -b_der; diff --git a/src/objects/object_general.hpp b/src/objects/object_general.hpp index f9a9d25..49a4ca8 100644 --- a/src/objects/object_general.hpp +++ b/src/objects/object_general.hpp @@ -8,17 +8,17 @@ class ObjectGeneral{ public: - BoundaryCondition bc; + BC bc; T friction; std::string name; - ObjectGeneral(BoundaryCondition bc, T friction, std::string name) : bc(bc), friction(friction), name(name) {} + ObjectGeneral(BC bc, T friction, std::string name) : bc(bc), friction(friction), name(name) {} virtual ~ObjectGeneral(){} - virtual bool inside(const TV& X_in) = 0; + virtual bool inside(const TV& X_in) const = 0; - virtual TV normal(const TV& X_in) = 0; + virtual TV normal(const TV& X_in) const = 0; }; diff --git a/src/objects/object_ground.hpp b/src/objects/object_ground.hpp index cb41730..55b8d10 100644 --- a/src/objects/object_ground.hpp +++ b/src/objects/object_ground.hpp @@ -11,9 +11,9 @@ class ObjectGround : public ObjectGeneral { ~ObjectGround(){} - ObjectGround(BoundaryCondition bc_in = NOSLIP, T friction_in = 0.0, std::string name_in = "", T y_ground_in = 0) : ObjectGeneral(bc_in, friction_in, name_in), y_ground(y_ground_in) {} + ObjectGround(BC bc_in = BC::NoSlip, T friction_in = 0.0, std::string name_in = "", T y_ground_in = 0) : ObjectGeneral(bc_in, friction_in, name_in), y_ground(y_ground_in) {} - bool inside(const TV& X_in) override { + bool inside(const TV& X_in) const override { T y = X_in(1); @@ -24,7 +24,7 @@ class ObjectGround : public ObjectGeneral { } - TV normal(const TV& X_in) override { + TV normal(const TV& X_in) const override { TV n; diff --git a/src/objects/object_plate.hpp b/src/objects/object_plate.hpp index 52cba36..a01b414 100644 --- a/src/objects/object_plate.hpp +++ b/src/objects/object_plate.hpp @@ -10,7 +10,7 @@ class ObjectPlate{ #ifdef THREEDIM - ObjectPlate(T pos_object, PlateType plate_type, BoundaryCondition bc = NOSLIP, T friction = 0, T pos_lower = -1e15, T pos_upper = 1e15, T vx_object = 0.0, T vy_object = 0.0, T vz_object = 0.0, T vmin_factor = 1.0, T load_factor = 0.0, std::string name = "") : + ObjectPlate(T pos_object, PlateType plate_type, BC bc = BC::NoSlip, T friction = 0, T pos_lower = -1e15, T pos_upper = 1e15, T vx_object = 0.0, T vy_object = 0.0, T vz_object = 0.0, T vmin_factor = 1.0, T load_factor = 0.0, std::string name = "") : pos_object(pos_object), plate_type(plate_type), bc(bc), @@ -27,43 +27,43 @@ class ObjectPlate{ load_factor(load_factor), name(name) {} - bool inside(const TV& X_in){ - if (plate_type == left){ + bool inside(const TV& X_in) const { + if (plate_type == PlateType::left){ if (X_in(1) < pos_upper && X_in(1) > pos_lower && (X_in(0) - pos_object) <= 0){ return true; } else{ return false; } } - else if (plate_type == right){ + else if (plate_type == PlateType::right){ if (X_in(1) < pos_upper && X_in(1) > pos_lower && (pos_object - X_in(0)) <= 0){ return true; } else{ return false; } } - else if (plate_type == bottom){ + else if (plate_type == PlateType::bottom){ if (X_in(0) < pos_upper && X_in(0) > pos_lower && (X_in(1) - pos_object) <= 0){ return true; } else{ return false; } } - else if (plate_type == top){ + else if (plate_type == PlateType::top){ if (X_in(0) < pos_upper && X_in(0) > pos_lower && (pos_object - X_in(1)) <= 0){ return true; } else{ return false; } } - else if (plate_type == back){ + else if (plate_type == PlateType::back){ if (X_in(2) - pos_object <= 0){ return true; } else{ return false; } } - else if (plate_type == front){ + else if (plate_type == PlateType::front){ if (pos_object - X_in(2) <= 0){ return true; } else{ @@ -80,7 +80,7 @@ class ObjectPlate{ #else // TWODIM - ObjectPlate(T pos_object, PlateType plate_type, BoundaryCondition bc = NOSLIP, T friction = 0, T pos_lower = -1e15, T pos_upper = 1e15, T vx_object = 0.0, T vy_object = 0.0, T vmin_factor = 1.0, T load_factor = 0.0, std::string name = "") : + ObjectPlate(T pos_object, PlateType plate_type, BC bc = BC::NoSlip, T friction = 0, T pos_lower = -1e15, T pos_upper = 1e15, T vx_object = 0.0, T vy_object = 0.0, T vmin_factor = 1.0, T load_factor = 0.0, std::string name = "") : pos_object(pos_object), plate_type(plate_type), bc(bc), @@ -95,29 +95,29 @@ class ObjectPlate{ load_factor(load_factor), name(name) {} - bool inside(const TV& X_in){ - if (plate_type == left){ + bool inside(const TV& X_in) const { + if (plate_type == PlateType::left){ if (X_in(1) < pos_upper && X_in(1) > pos_lower && (X_in(0) - pos_object) <= 0){ return true; } else{ return false; } } - else if (plate_type == right){ + else if (plate_type == PlateType::right){ if (X_in(1) < pos_upper && X_in(1) > pos_lower && (pos_object - X_in(0)) <= 0){ return true; } else{ return false; } } - else if (plate_type == bottom){ + else if (plate_type == PlateType::bottom){ if (X_in(0) < pos_upper && X_in(0) > pos_lower && (X_in(1) - pos_object) <= 0){ return true; } else{ return false; } } - else if (plate_type == top){ + else if (plate_type == PlateType::top){ if (X_in(0) < pos_upper && X_in(0) > pos_lower && (pos_object - X_in(1)) <= 0){ return true; } else{ @@ -135,34 +135,25 @@ class ObjectPlate{ void move(T dt, T frame_dt, T time){ - T load_time = load_factor * frame_dt; - if (time < load_time) { - vx_object = vx_object_original / vmin_factor; - vy_object = vy_object_original / vmin_factor; - #ifdef THREEDIM - vz_object = vz_object_original / vmin_factor; - #endif - } - else{ - vx_object = vx_object_original; - vy_object = vy_object_original; - #ifdef THREEDIM - vz_object = vz_object_original; - #endif - } + T factor = (time < load_factor * frame_dt) ? (1.0 / vmin_factor) : 1.0; + vx_object = vx_object_original * factor; + vy_object = vy_object_original * factor; + #ifdef THREEDIM + vz_object = vz_object_original * factor; + #endif - if (plate_type == left || plate_type == right){ + if (plate_type == PlateType::left || plate_type == PlateType::right){ pos_object += dt * vx_object; pos_upper += dt * vy_object; pos_lower += dt * vy_object; } - else if (plate_type == bottom || plate_type == top){ + else if (plate_type == PlateType::bottom || plate_type == PlateType::top){ pos_object += dt * vy_object; pos_upper += dt * vx_object; pos_lower += dt * vx_object; } #ifdef THREEDIM - else if (plate_type == back || plate_type == front){ + else if (plate_type == PlateType::back || plate_type == PlateType::front){ pos_object += dt * vz_object; } #endif @@ -173,7 +164,7 @@ class ObjectPlate{ T pos_object; PlateType plate_type; - BoundaryCondition bc; + BC bc; T friction; T pos_upper; diff --git a/src/objects/object_ramp.hpp b/src/objects/object_ramp.hpp index 8798732..4e99ab7 100644 --- a/src/objects/object_ramp.hpp +++ b/src/objects/object_ramp.hpp @@ -10,9 +10,9 @@ class ObjectRamp : public ObjectGeneral { ~ObjectRamp(){} - ObjectRamp(BoundaryCondition bc_in, T friction_in, std::string name_in = "") : ObjectGeneral(bc_in, friction_in, name_in) {} + ObjectRamp(BC bc_in, T friction_in, std::string name_in = "") : ObjectGeneral(bc_in, friction_in, name_in) {} - bool inside(const TV& X_in) override { + bool inside(const TV& X_in) const override { T x = X_in(0); T y = X_in(1); @@ -25,10 +25,10 @@ class ObjectRamp : public ObjectGeneral { } - TV normal(const TV& X_in) override { + TV normal(const TV& X_in) const override { T x = X_in(0); - TV n; + TV n = TV::Zero(); T tmp = std::cosh(-100*x); T b_der = 0.1*(-100)/(tmp*tmp); diff --git a/src/objects/object_silo.hpp b/src/objects/object_silo.hpp index 0ea37c8..ca7ac3f 100644 --- a/src/objects/object_silo.hpp +++ b/src/objects/object_silo.hpp @@ -11,9 +11,9 @@ class ObjectSilo : public ObjectGeneral { ~ObjectSilo(){} - ObjectSilo(BoundaryCondition bc_in, T friction_in, std::string name_in = "", T cut_in = -1) : ObjectGeneral(bc_in, friction_in, name_in), cut(cut_in) {} + ObjectSilo(BC bc_in, T friction_in, std::string name_in = "", T cut_in = -1) : ObjectGeneral(bc_in, friction_in, name_in), cut(cut_in) {} - bool inside(const TV& X_in) override { + bool inside(const TV& X_in) const override { T x = X_in(0); T y = X_in(1); @@ -34,7 +34,7 @@ class ObjectSilo : public ObjectGeneral { } // end inside - TV normal(const TV& X_in) override { + TV normal(const TV& X_in) const override { TV n; diff --git a/src/objects/object_vdb.hpp b/src/objects/object_vdb.hpp index dfb630d..ccd4876 100644 --- a/src/objects/object_vdb.hpp +++ b/src/objects/object_vdb.hpp @@ -26,7 +26,7 @@ class ObjectVdb : public ObjectGeneral { ~ObjectVdb(){} - ObjectVdb(std::string filename, BoundaryCondition bc_in = NOSLIP, T friction_in = 0.0, std::string name_in = "") : ObjectGeneral(bc_in, friction_in, name_in) { + ObjectVdb(std::string filename, BC bc_in = BC::NoSlip, T friction_in = 0.0, std::string name_in = "") : ObjectGeneral(bc_in, friction_in, name_in) { openvdb::io::File file(filename); file.open(); @@ -42,7 +42,7 @@ class ObjectVdb : public ObjectGeneral { grad_phi = mg.process(); } - bool inside(const TV& X_in) override { + bool inside(const TV& X_in) const override { int dim = X_in.size(); Eigen::Matrix X; X.setZero(); @@ -56,7 +56,7 @@ class ObjectVdb : public ObjectGeneral { return ((T)phi <= 0); } - TV normal(const TV& X_in) override { + TV normal(const TV& X_in) const override { int dim = X_in.size(); Eigen::Matrix X; X.setZero(); @@ -76,7 +76,7 @@ class ObjectVdb : public ObjectGeneral { return TV::Zero(); } - void bounds(TV& min_bbox, TV& max_bbox){ + void bounds(TV& min_bbox, TV& max_bbox) const { int dim = min_bbox.size(); min_bbox.setZero(); max_bbox.setZero(); diff --git a/src/plasticity_helpers/mcc_rma_explicit.hpp b/src/plasticity_helpers/mcc_rma_explicit.hpp index 4c42581..791d07a 100644 --- a/src/plasticity_helpers/mcc_rma_explicit.hpp +++ b/src/plasticity_helpers/mcc_rma_explicit.hpp @@ -5,7 +5,7 @@ #include "../tools.hpp" -bool MCCRMAExplicit(T& p, T& q, int& exit, T M, T p0, T beta, T mu, T K, T rma_prefac) +bool MCCRMAExplicit(T& p, T& q, int& exit, T M, T p0, T beta, T mu, T K, T f_mu_prefac) { // T y = M * M * (p - p0) * (p + beta * p0) + (1 + 2 * beta) * (q * q); T y = M * M * (p - p0) * (p + beta * p0) + (q * q); @@ -30,8 +30,8 @@ bool MCCRMAExplicit(T& p, T& q, int& exit, T M, T p0, T beta, T mu, T K, T rma_p // T l = 2*q * (2*beta+1); T l = 2*q; - T r1 = pt - p - K * delta_gamma * k; - T r2 = qt - q - rma_prefac*mu * delta_gamma * l; + T r1 = pt - p - K * delta_gamma * k; + T r2 = qt - q - f_mu_prefac * delta_gamma * l; if ( iter > 4 && std::abs(y) < 1e-3 && std::abs(r1) < 1e-3 && std::abs(r2) < 1e-3 ){ break; @@ -57,11 +57,11 @@ bool MCCRMAExplicit(T& p, T& q, int& exit, T M, T p0, T beta, T mu, T K, T rma_p } } - T J11 = -1 - K *delta_gamma*dkdp; + T J11 = -1 - K *delta_gamma*dkdp; T J13 = -K*k; - T J22 = -1 - rma_prefac*mu*delta_gamma*dldq; - T J23 = -rma_prefac*mu*l; + T J22 = -1 - f_mu_prefac*delta_gamma*dldq; + T J23 = -f_mu_prefac*l; T J31 = k; T J32 = l; diff --git a/src/plasticity_helpers/mcc_rma_explicit_onevar.hpp b/src/plasticity_helpers/mcc_rma_explicit_onevar.hpp index 0eea490..7133e3f 100644 --- a/src/plasticity_helpers/mcc_rma_explicit_onevar.hpp +++ b/src/plasticity_helpers/mcc_rma_explicit_onevar.hpp @@ -5,7 +5,7 @@ #include "../tools.hpp" -bool MCCRMAExplicitOnevar(T& p, T& q, int& exit, T M, T p0, T beta, T mu, T K, T rma_prefac) +bool MCCRMAExplicitOnevar(T& p, T& q, int& exit, T M, T p0, T beta, T mu, T K) /** **** WRITTEN BY TOBIAS VERHEIJEN **** * @@ -17,7 +17,6 @@ bool MCCRMAExplicitOnevar(T& p, T& q, int& exit, T M, T p0, T beta, T mu, T K, T * @param beta: Cohesion parameter * @param mu: G = Lame parameter * @param K: Bulk modulus, (lambda + 2*mu/dim) - * @param rma_prefac: precomputed value, either 2*sqrt(3) or 1. depending on using von_mises_q or not * @return true if deformation is plastic, false if deformation is not plastic */ { diff --git a/src/plasticity_helpers/mcc_rma_implicit_exponential.hpp b/src/plasticity_helpers/mcc_rma_implicit_exponential.hpp index 00dd67f..2af8c82 100644 --- a/src/plasticity_helpers/mcc_rma_implicit_exponential.hpp +++ b/src/plasticity_helpers/mcc_rma_implicit_exponential.hpp @@ -5,7 +5,7 @@ #include "../tools.hpp" -bool MCCRMAImplicitExponential(T& p, T& q, int& exit, T M, T p00, T beta, T mu, T K, T xi, T rma_prefac, T epv) +bool MCCRMAImplicitExponential(T& p, T& q, int& exit, T M, T p00, T beta, T mu, T K, T xi, T f_mu_prefac, T epv) { T p0 = std::max(T(1e-2), p00 * std::exp(-xi*epv)); // T y = M * M * (p - p0) * (p + beta * p0) + (1 + 2 * beta) * (q * q); @@ -34,8 +34,8 @@ bool MCCRMAImplicitExponential(T& p, T& q, int& exit, T M, T p00, T beta, T mu, // T dydq = 2*q * (2*beta+1); T dydq = 2*q; - T r1 = pt - p - K * delta_gamma * dydp; - T r2 = qt - q - rma_prefac*mu * delta_gamma * dydq; + T r1 = pt - p - K * delta_gamma * dydp; + T r2 = qt - q - f_mu_prefac * delta_gamma * dydq; if ( iter > 4 && std::abs(y) < 1e-3 && std::abs(r1) < 1e-3 && std::abs(r2) < 1e-3 ){ break; @@ -60,11 +60,11 @@ bool MCCRMAImplicitExponential(T& p, T& q, int& exit, T M, T p00, T beta, T mu, } } - T J11 = -1 - K *delta_gamma*ddydpp; + T J11 = -1 - K *delta_gamma*ddydpp; T J13 = -K*dydp; - T J22 = -1 - rma_prefac*mu*delta_gamma*ddydqq; - T J23 = -rma_prefac*mu*dydq; + T J22 = -1 - f_mu_prefac*delta_gamma*ddydqq; + T J23 = -f_mu_prefac*dydq; T J31 = dydp; T J32 = dydq; diff --git a/src/plasticity_helpers/mcc_rma_implicit_exponential_onevar.hpp b/src/plasticity_helpers/mcc_rma_implicit_exponential_onevar.hpp index 0b7a1ce..46c948f 100644 --- a/src/plasticity_helpers/mcc_rma_implicit_exponential_onevar.hpp +++ b/src/plasticity_helpers/mcc_rma_implicit_exponential_onevar.hpp @@ -6,7 +6,7 @@ #include "../tools.hpp" #include "limited_search_exponential.hpp" -bool MCCRMAImplicitExponentialOnevar(T& p, T& q, int& exit, T M, T p00, T beta, T mu, T K, T xi, T rma_prefac, T epv) +bool MCCRMAImplicitExponentialOnevar(T& p, T& q, int& exit, T M, T p00, T beta, T mu, T K, T xi, T epv) /** **** WRITTEN BY TOBIAS VERHEIJEN **** * @@ -30,7 +30,6 @@ bool MCCRMAImplicitExponentialOnevar(T& p, T& q, int& exit, T M, T p00, T beta, * @param K: Bulk modulus, (lambda + 2*mu/dim) * @param xi: xi value for simulation * @param epv: Plastic Volumetric Hencky Strain - * @param rma_prefac: precomputed value, either 2*sqrt(3) or 1. depending on using von_mises_q or not * @return true if deformation is plastic, false if deformation is not plastic */ { @@ -45,9 +44,15 @@ bool MCCRMAImplicitExponentialOnevar(T& p, T& q, int& exit, T M, T p00, T beta, return false; } - // store initial values and rescale - T scale_factor = 1./p0_t; - T p0 = 1.; p00 *= scale_factor; p *= scale_factor; q *= scale_factor; mu *= scale_factor; K *= scale_factor; + // ALT 1: Rescale + T scale_factor = 1./p0_t; + T p0 = 1.; + + // ALT 2: Do not rescale + // T scale_factor = 1; + // T p0 = p0_t; + + p00 *= scale_factor; p *= scale_factor; q *= scale_factor; mu *= scale_factor; K *= scale_factor; T p_t = p; T q_t = q; T epv_t = epv; @@ -87,33 +92,27 @@ bool MCCRMAImplicitExponentialOnevar(T& p, T& q, int& exit, T M, T p00, T beta, y = Msq * (p - p0) * (p + beta * p0) + damping_factor * q * q; // Check for convergence, iteration > 3 check is done to ensure some displacement is measured. - if (iter > 3 && abs(y) < 1.e-3 * scale_factor){ + if (iter > 5 && abs(y) < 1.e-4){ break; } if (iter == max_iter - 1) { // did not break loop if (p0 > 1.e-3 * scale_factor) { debug("RMA: FATAL did not exit loop at iteration = ", iter, ", iteration = ", iter); - debug(iter, ": r1 = ", y); - debug(iter, ": y = ", y); - debug(iter, ": p0 = ", p0); - debug(iter, ": p0_t = ", p0_t); - debug(iter, ": K = ", K); - debug(iter, ": G = ", mu); - debug(iter, ": check value = ", 1.e-3 * scale_factor); - debug(iter, ": pt = ", p_t / scale_factor); - debug(iter, ": qt = ", q_t / scale_factor); - debug(iter, ": scale factor = ", scale_factor); - debug(iter, ": epvt = ", epv_t); - debug(iter, ": epv = ", epv); - debug(iter, ": p = ", p / scale_factor); - debug(iter, ": q = ", q / scale_factor); - - exit = 1; - } else { // p0 too small - p = 1e-15; - q = 1e-15; - break; + debug("Consider changing the threshold for breaking the N-R loop"); + debug("y = ", y); + debug("p0 = ", p0); + debug("K = ", K); + debug("G = ", mu); + debug("pt = ", p_t); + debug("qt = ", q_t); + debug("epvt = ", epv_t); + debug("epv = ", epv); + debug("p = ", p); + debug("q = ", q); + debug("p0_t = ", p0_t); + debug("scale factor = ", scale_factor); + // exit = 1; } } @@ -133,20 +132,19 @@ bool MCCRMAImplicitExponentialOnevar(T& p, T& q, int& exit, T M, T p00, T beta, epv = limited_search_exponential(epv, epv_t, y / dy_depv, xi, K, p_t, p00, beta, right); if (epv == std::numeric_limits::infinity()) { - debug("RMA: FATAL no valid epv value found = ", iter, ", iteration = ", iter); - debug(iter, ": r1 = ", y); - debug(iter, ": y = ", y); - debug(iter, ": p0 = ", p0); - debug(iter, ": p0_t = ", p0_t); - debug(iter, ": pt = ", p_t); - debug(iter, ": qt = ", q_t); - debug(iter, ": epvt = ", epv_t); - debug(iter, ": epv = ", epv); - debug(iter, ": p = ", p); - debug(iter, ": q = ", q); - debug(iter, ": K = ", K); - debug(iter, ": mu = ", mu); - + debug("RMA: FATAL no valid epv value found = ", iter); + debug("y = ", y); + debug("p0 = ", p0); + debug("K = ", K); + debug("G = ", mu); + debug("pt = ", p_t); + debug("qt = ", q_t); + debug("epvt = ", epv_t); + debug("epv = ", epv); + debug("p = ", p); + debug("q = ", q); + debug("p0_t = ", p0_t); + debug("scale factor = ", scale_factor); exit = 1; break; } diff --git a/src/plasticity_helpers/mcc_rma_implicit_sinh_onevar.hpp b/src/plasticity_helpers/mcc_rma_implicit_sinh_onevar.hpp index f69950b..3a3fea3 100644 --- a/src/plasticity_helpers/mcc_rma_implicit_sinh_onevar.hpp +++ b/src/plasticity_helpers/mcc_rma_implicit_sinh_onevar.hpp @@ -6,7 +6,7 @@ #include "../tools.hpp" #include "limited_search_sinh.hpp" -bool MCCRMAImplicitSinhOnevar(T& p, T& q, int& exit, T M, T p00, T beta, T mu, T K, T xi, T rma_prefac, T epv) +bool MCCRMAImplicitSinhOnevar(T& p, T& q, int& exit, T M, T p00, T beta, T mu, T K, T xi, T epv) /** **** WRITTEN BY TOBIAS VERHEIJEN **** * @@ -30,7 +30,6 @@ bool MCCRMAImplicitSinhOnevar(T& p, T& q, int& exit, T M, T p00, T beta, T mu, T * @param K: Bulk modulus, (lambda + 2*mu/dim) * @param xi: xi value for simulation * @param epv: Plastic Volumetric Hencky Strain - * @param rma_prefac: precomputed value, either 2*sqrt(3) or 1. depending on using von_mises_q or not respectively * @return true if deformation is plastic, false if deformation is not plastic */ { @@ -90,31 +89,26 @@ bool MCCRMAImplicitSinhOnevar(T& p, T& q, int& exit, T M, T p00, T beta, T mu, T y = Msq * (p - p0) * (p + beta * p0) + damping_factor * q * q; // Check for convergence, iteration > 3 check is done to ensure some displacement is measured. - if (iter > 3 && std::abs(y) < 1e-3 * scale_factor) { + if (iter > 5 && std::abs(y) < 1e-4) { break; } if (iter == max_iter - 1) { // did not break loop if (p0 > p0_min) { debug("RMA: FATAL did not exit loop at iteration = ", iter, ", iteration = ", iter); - debug(iter, ": r1 = ", y); - debug(iter, ": y = ", y); - debug(iter, ": p0 = ", p0); - debug(iter, ": p0_t = ", p0_t); - debug(iter, ": pt = ", p_t); - debug(iter, ": qt = ", q_t); - debug(iter, ": epv_t = ", epv_t); - debug(iter, ": epv = ", epv); - debug(iter, ": p = ", p); - debug(iter, ": q = ", q); - debug(iter, ": K = ", K); - debug(iter, ": mu = ", mu); - + debug("y = ", y); + debug("p0 = ", p0); + debug("K = ", K); + debug("G = ", mu); + debug("pt = ", p_t); + debug("qt = ", q_t); + debug("epvt = ", epv_t); + debug("epv = ", epv); + debug("p = ", p); + debug("q = ", q); + debug("p0_t = ", p0_t); + debug("scale factor = ", scale_factor); exit = 1; - } else { // p0 too small - p = 1e-15; - q = 1e-15; - break; } } @@ -156,7 +150,6 @@ bool MCCRMAImplicitSinhOnevar(T& p, T& q, int& exit, T M, T p00, T beta, T mu, T debug(iter, ": q = ", q); debug(iter, ": K = ", K); debug(iter, ": mu = ", mu); - exit = 1; break; } diff --git a/src/sampling/sampling_particles.hpp b/src/sampling/sampling_particles.hpp index dc8cf93..a2bdb48 100644 --- a/src/sampling/sampling_particles.hpp +++ b/src/sampling/sampling_particles.hpp @@ -122,7 +122,7 @@ } // end sampleParticles -#endif // DIMENSION +#endif #endif // SAMPLING_PARTICLES_HPP diff --git a/src/sampling/sampling_particles_vdb.hpp b/src/sampling/sampling_particles_vdb.hpp index 2723b91..37d089b 100644 --- a/src/sampling/sampling_particles_vdb.hpp +++ b/src/sampling/sampling_particles_vdb.hpp @@ -14,7 +14,7 @@ template void sampleParticlesFromVdb(S& sim, ObjectVdb& obj, T kRadius, T ppc = 8) #else // TWODIM void sampleParticlesFromVdb(S& sim, ObjectVdb& obj, T kRadius, T ppc = 6) -#endif // DIMENSION +#endif { debug("Sampling particles from VDB..."); @@ -38,6 +38,9 @@ void sampleParticlesFromVdb(S& sim, ObjectVdb& obj, T kRadius, T ppc = 6) sim.particle_volume = sim.dx * sim.dx * sim.dx / ppc; sim.particle_mass = sim.rho * sim.particle_volume; + sim.Lx = L(0); + sim.Ly = L(1); + sim.Lz = L(2); #else // TWODIM debug(" Min corner: ", min_corner(0), ", ", min_corner(1)); debug(" Max corner: ", max_corner(0), ", ", max_corner(1)); @@ -49,7 +52,10 @@ void sampleParticlesFromVdb(S& sim, ObjectVdb& obj, T kRadius, T ppc = 6) sim.dx = std::sqrt(ppc / T(square_samples.size()) * L(0)*L(1)); sim.particle_volume = sim.dx * sim.dx / ppc; sim.particle_mass = sim.rho * sim.particle_volume; - #endif // DIMENSION + + sim.Lx = L(0); + sim.Ly = L(1); + #endif debug(" Number of square samples: ", square_samples.size()); debug(" dx set to ", sim.dx); @@ -61,7 +67,7 @@ void sampleParticlesFromVdb(S& sim, ObjectVdb& obj, T kRadius, T ppc = 6) TV point(square_samples[p][0], square_samples[p][1], square_samples[p][2]); #else // TWODIM TV point(square_samples[p][0], square_samples[p][1]); - #endif // DIMENSION + #endif if ( obj.inside(point) ){ samples.push_back(point); diff --git a/src/simulation/boundary_collision.cpp b/src/simulation/boundary_collision.cpp index 0ac596c..227815e 100644 --- a/src/simulation/boundary_collision.cpp +++ b/src/simulation/boundary_collision.cpp @@ -10,16 +10,16 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // New positions Xi += dt * vi; // recall Xi was passed by value - for(auto obj : objects) { + for(auto& obj : objects) { bool colliding = obj->inside(Xi); if (colliding) { TV v_rel = vi_orig; - if (obj->bc == NOSLIP) { + if (obj->bc == BC::NoSlip) { v_rel.setZero(); - } // end NOSLIP + } // end BC::NoSlip - else if (obj->bc == SLIPFREE) { + else if (obj->bc == BC::SlipFree) { TV n = obj->normal(Xi); T dot = v_rel.dot(n); if (dot < 0){ // if moving towards object @@ -39,9 +39,9 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ } // end non-zero friction } // end moving towards object - } // end SLIPFREE + } // end BC::SlipFree - else if (obj->bc == SLIPSTICK) { + else if (obj->bc == BC::SlipStick) { TV n = obj->normal(Xi); T dot = v_rel.dot(n); @@ -61,7 +61,7 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ } // end non-zero friction } // end moving towards object - } // end SLIPSTICK + } // end BC::SlipStick else { debug("INVALID BOUNDARY CONDITION!!!"); @@ -81,26 +81,26 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ #ifdef THREEDIM - for (ObjectPlate &obj : plates) { - bool colliding = obj.inside(Xi); + for (auto& obj : plates) { + bool colliding = obj->inside(Xi); if (colliding) { - T vx_rel = vi_orig(0) - obj.vx_object; - T vy_rel = vi_orig(1) - obj.vy_object; - T vz_rel = vi_orig(2) - obj.vz_object; + T vx_rel = vi_orig(0) - obj->vx_object; + T vy_rel = vi_orig(1) - obj->vy_object; + T vz_rel = vi_orig(2) - obj->vz_object; - if (obj.bc == NOSLIP) { + if (obj->bc == BC::NoSlip) { vx_rel = 0; vy_rel = 0; vz_rel = 0; - } // end NOSLIP + } // end BC::NoSlip - else if (obj.bc == SLIPSTICK) { + else if (obj->bc == BC::SlipStick) { - T friction = obj.friction; + T friction = obj->friction; if (use_mibf) friction = grid.friction[index]; - if (obj.plate_type == top || obj.plate_type == bottom){ + if (obj->plate_type == PlateType::top || obj->plate_type == PlateType::bottom){ // tangential velocity is the (x,z) components T vel_t = std::sqrt(vx_rel*vx_rel + vz_rel*vz_rel); @@ -121,8 +121,8 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // normal component (y) must be set to zero vy_rel = 0; - } // end top and bottom plate - else if (obj.plate_type == left || obj.plate_type == right){ + } // end PlateType::top and PlateType::bottom plate + else if (obj->plate_type == PlateType::left || obj->plate_type == PlateType::right){ // tangential velocity is the (y,z) components T vel_t = std::sqrt(vy_rel*vy_rel + vz_rel*vz_rel); @@ -143,8 +143,8 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // normal component (x) must be set to zero vx_rel = 0; - } // end left or right plate - else if (obj.plate_type == front || obj.plate_type == back){ + } // end PlateType::left or PlateType::right plate + else if (obj->plate_type == PlateType::front || obj->plate_type == PlateType::back){ // tangential velocity is the (x,y) components T vel_t = std::sqrt(vx_rel*vx_rel + vy_rel*vy_rel); @@ -165,18 +165,18 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // normal component (z) must be set to zero vz_rel = 0; - } // end front or back plate - } // end SLIPSTICK + } // end PlateType::front or PlateType::back plate + } // end BC::SlipStick - else if (obj.bc == SLIPFREE) { + else if (obj->bc == BC::SlipFree) { - T friction = obj.friction; + T friction = obj->friction; if (use_mibf) friction = grid.friction[index]; - if ((obj.plate_type == top && vy_rel > 0) || (obj.plate_type == bottom && vy_rel < 0)){ + if ((obj->plate_type == PlateType::top && vy_rel > 0) || (obj->plate_type == PlateType::bottom && vy_rel < 0)){ // tangential velocity is the (x,z) components T vel_t = std::sqrt(vx_rel*vx_rel + vz_rel*vz_rel); T fric_vel_n = friction * std::abs(vy_rel); @@ -196,8 +196,8 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // normal component (y) must be set to zero vy_rel = 0; - } // end top and bottom plate - else if ((obj.plate_type == left && vx_rel < 0) || (obj.plate_type == right && vx_rel > 0)){ + } // end PlateType::top and PlateType::bottom plate + else if ((obj->plate_type == PlateType::left && vx_rel < 0) || (obj->plate_type == PlateType::right && vx_rel > 0)){ // tangential velocity is the (y,z) components T vel_t = std::sqrt(vy_rel*vy_rel + vz_rel*vz_rel); T fric_vel_n = friction * std::abs(vx_rel); @@ -217,8 +217,8 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // normal component (x) must be set to zero vx_rel = 0; - } // end left or right plate - else if ((obj.plate_type == front && vz_rel > 0) || (obj.plate_type == back && vz_rel < 0 )){ + } // end PlateType::left or PlateType::right plate + else if ((obj->plate_type == PlateType::front && vz_rel > 0) || (obj->plate_type == PlateType::back && vz_rel < 0 )){ // tangential velocity is the (x,y) components T vel_t = std::sqrt(vx_rel*vx_rel + vy_rel*vy_rel); T fric_vel_n = friction * std::abs(vz_rel); @@ -238,8 +238,8 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // normal component (z) must be set to zero vz_rel = 0; - } // end front or back plate - } // end SLIPFREE + } // end PlateType::front or PlateType::back plate + } // end BC::SlipFree @@ -249,9 +249,9 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ return; } - vi(0) = vx_rel + obj.vx_object; - vi(1) = vy_rel + obj.vy_object; - vi(2) = vz_rel + obj.vz_object; + vi(0) = vx_rel + obj->vx_object; + vi(1) = vy_rel + obj->vy_object; + vi(2) = vz_rel + obj->vz_object; // update velocity copy before next iteration vi_orig = vi; // Comment this line to enforce ordering of objects (i.e., use only last object in list) @@ -264,24 +264,24 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ #else // TWODIM - for (ObjectPlate &obj : plates) { - bool colliding = obj.inside(Xi); + for (auto& obj : plates) { + bool colliding = obj->inside(Xi); if (colliding) { - T vx_rel = vi_orig(0) - obj.vx_object; - T vy_rel = vi_orig(1) - obj.vy_object; + T vx_rel = vi_orig(0) - obj->vx_object; + T vy_rel = vi_orig(1) - obj->vy_object; - if (obj.bc == NOSLIP) { + if (obj->bc == BC::NoSlip) { vx_rel = 0; vy_rel = 0; - } // end NOSLIP + } // end BC::NoSlip - else if (obj.bc == SLIPSTICK) { + else if (obj->bc == BC::SlipStick) { - T friction = obj.friction; + T friction = obj->friction; if (use_mibf) friction = grid.friction[index]; - if (obj.plate_type == top || obj.plate_type == bottom){ + if (obj->plate_type == PlateType::top || obj->plate_type == PlateType::bottom){ // tangential velocity is the (x) component and must be changed if (friction == 0){ // Do nothing @@ -299,8 +299,8 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // normal component (y) must be set to zero vy_rel = 0; - } // end top and bottom plate - else if (obj.plate_type == left || obj.plate_type == right){ + } // end PlateType::top and PlateType::bottom plate + else if (obj->plate_type == PlateType::left || obj->plate_type == PlateType::right){ // tangential velocity is the (y) component and must be changed if (friction == 0){ @@ -319,19 +319,19 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // normal component (x) must be set to zero vx_rel = 0; - } // end left or right plate + } // end PlateType::left or PlateType::right plate - } // end SLIPSTICK + } // end BC::SlipStick - else if (obj.bc == SLIPFREE) { + else if (obj->bc == BC::SlipFree) { - T friction = obj.friction; + T friction = obj->friction; if (use_mibf) friction = grid.friction[index]; - if ((obj.plate_type == top && vy_rel > 0) || (obj.plate_type == bottom && vy_rel < 0)){ + if ((obj->plate_type == PlateType::top && vy_rel > 0) || (obj->plate_type == PlateType::bottom && vy_rel < 0)){ // tangential velocity is the (x) component and must be changed if (friction == 0){ // Do nothing @@ -349,8 +349,8 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // normal component (y) must be set to zero vy_rel = 0; - } // end top and bottom plate - else if ((obj.plate_type == left && vx_rel < 0) || (obj.plate_type == right && vx_rel > 0)){ + } // end PlateType::top and PlateType::bottom plate + else if ((obj->plate_type == PlateType::left && vx_rel < 0) || (obj->plate_type == PlateType::right && vx_rel > 0)){ // tangential velocity is the (y) component and must be changed if (friction == 0){ // Do nothing @@ -368,8 +368,8 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ // normal component (x) must be set to zero vx_rel = 0; - } // end left or right plate - } // end SLIPFREE + } // end PlateType::left or PlateType::right plate + } // end BC::SlipFree @@ -379,8 +379,8 @@ void Simulation::boundaryCollision(int index, TV Xi, TV& vi){ return; } - vi(0) = vx_rel + obj.vx_object; - vi(1) = vy_rel + obj.vy_object; + vi(0) = vx_rel + obj->vx_object; + vi(1) = vy_rel + obj->vy_object; // update velocity copy before next iteration vi_orig = vi; // Comment this line to enforce ordering of objects (i.e., use only last object in list) diff --git a/src/simulation/explicit_euler_update.cpp b/src/simulation/explicit_euler_update.cpp index 5533107..d7b003a 100644 --- a/src/simulation/explicit_euler_update.cpp +++ b/src/simulation/explicit_euler_update.cpp @@ -16,16 +16,16 @@ void Simulation::explicitEulerUpdate(){ { std::vector grid_force_local(grid_nodes, TV::Zero()); - #pragma omp for + #pragma omp for nowait for(int p = 0; p < Np; p++){ TM Fe = particles.F[p]; TM dPsidF; - if (elastic_model == NeoHookean){ + if (elastic_model == ElasticModel::NeoHookean){ dPsidF = NeoHookeanPiola(Fe); } - else if (elastic_model == Hencky){ // St Venant Kirchhoff with Hencky strain + else if (elastic_model == ElasticModel::Hencky){ // St Venant Kirchhoff with Hencky strain dPsidF = HenckyPiola(Fe); } else{ @@ -33,7 +33,6 @@ void Simulation::explicitEulerUpdate(){ } TM tau = dPsidF * Fe.transpose(); - particles.tau[p] = tau; TV xp = particles.x[p]; unsigned int i_base = std::max(0, int(std::floor((xp(0)-grid.xc)*one_over_dx)) - 1); // i_base = std::min(i_base, Nx-4); // the subtraction of one is valid for both quadratic and cubic splines diff --git a/src/simulation/g2p.cpp b/src/simulation/g2p.cpp index b00469a..7686313 100644 --- a/src/simulation/g2p.cpp +++ b/src/simulation/g2p.cpp @@ -20,7 +20,7 @@ void Simulation::G2P(){ std::vector particles_flip_local(Np); std::fill( particles_flip_local.begin(), particles_flip_local.end(), TV::Zero() ); std::vector particles_Bmat_local(Np); std::fill( particles_Bmat_local.begin(), particles_Bmat_local.end(), TM::Zero() ); - #pragma omp for + #pragma omp for nowait for(int p = 0; p < Np; p++){ TV xp = particles.x[p]; TV vp = TV::Zero(); diff --git a/src/simulation/musl.cpp b/src/simulation/musl.cpp index 1450c70..585a7b9 100644 --- a/src/simulation/musl.cpp +++ b/src/simulation/musl.cpp @@ -10,7 +10,7 @@ void Simulation::MUSL(){ #endif - #pragma omp parallel for num_threads(n_threads) + #pragma omp parallel for schedule(static) num_threads(n_threads) for(int p=0; p grid_v_local(grid_nodes, TV::Zero() ); - #pragma omp for + #pragma omp for nowait for(int p = 0; p < Np; p++){ TV xp = particles.x[p]; unsigned int i_base = std::max(0, int(std::floor((xp(0)-grid.xc)*one_over_dx)) - 1); // i_base = std::min(i_base, Nx-4); // the subtraction of one is valid for both quadratic and cubic splines diff --git a/src/simulation/p2g.cpp b/src/simulation/p2g.cpp index 2189ff9..211b5f6 100644 --- a/src/simulation/p2g.cpp +++ b/src/simulation/p2g.cpp @@ -15,7 +15,7 @@ void Simulation::P2G(){ std::vector grid_mass_local(grid_nodes); std::vector grid_friction_local(grid_nodes); - #pragma omp for + #pragma omp for nowait for(int p = 0; p < Np; p++){ TV xp = particles.x[p]; unsigned int i_base = std::max(0, int(std::floor((xp(0)-grid.xc)*one_over_dx)) - 1); // i_base = std::min(i_base, Nx-4); // the subtraction of one is valid for both quadratic and cubic splines diff --git a/src/simulation/pbc.cpp b/src/simulation/pbc.cpp index d66aa2f..b49ff80 100644 --- a/src/simulation/pbc.cpp +++ b/src/simulation/pbc.cpp @@ -15,14 +15,12 @@ void Simulation::PBCAddParticles1D(){ particles.x.push_back(part_x+incr); particles.v.push_back(particles.v[p]); particles.F.push_back(particles.F[p]); - particles.tau.push_back(particles.tau[p]); particles.Bmat.push_back(particles.Bmat[p]); num_add_pbc_particles++; particles.x.push_back(part_x-incr); particles.v.push_back(particles.v[p]); particles.F.push_back(particles.F[p]); - particles.tau.push_back(particles.tau[p]); particles.Bmat.push_back(particles.Bmat[p]); num_add_pbc_particles++; @@ -54,7 +52,6 @@ void Simulation::PBCAddParticles(unsigned int safety_factor){ particles.x.push_back(part_x); particles.v.push_back(particles.v[p]); particles.F.push_back(particles.F[p]); - particles.tau.push_back(particles.tau[p]); particles.Bmat.push_back(particles.Bmat[p]); } num_add_pbc_particles++; @@ -69,7 +66,6 @@ void Simulation::PBCAddParticles(unsigned int safety_factor){ particles.x.push_back(part_x); particles.v.push_back(particles.v[p]); particles.F.push_back(particles.F[p]); - particles.tau.push_back(particles.tau[p]); particles.Bmat.push_back(particles.Bmat[p]); } num_add_pbc_particles++; @@ -95,7 +91,6 @@ void Simulation::PBCDelParticles(){ particles.x.erase( particles.x.end() -num_add_pbc_particles, particles.x.end() ); particles.v.erase( particles.v.end() -num_add_pbc_particles, particles.v.end() ); particles.F.erase( particles.F.end() -num_add_pbc_particles, particles.F.end() ); - particles.tau.erase( particles.tau.end() -num_add_pbc_particles, particles.tau.end() ); particles.Bmat.erase(particles.Bmat.end()-num_add_pbc_particles, particles.Bmat.end()); Np = particles.x.size(); diff --git a/src/simulation/plasticity.cpp b/src/simulation/plasticity.cpp index bb14a0a..b207156 100644 --- a/src/simulation/plasticity.cpp +++ b/src/simulation/plasticity.cpp @@ -9,11 +9,11 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & Fe_trial){ - if (plastic_model == NoPlasticity){ + if (plastic_model == PlasticModel::NoPlasticity){ // Do nothing } - else if (plastic_model == VM || plastic_model == DP || plastic_model == DPSoft || plastic_model == MCC || plastic_model == VMVisc || plastic_model == DPVisc || plastic_model == MCCVisc || plastic_model == DPMui || plastic_model == MCCMui){ + else if (plastic_model == PlasticModel::VM || plastic_model == PlasticModel::DP || plastic_model == PlasticModel::DPSoft || plastic_model == PlasticModel::MCC || plastic_model == PlasticModel::VMVisc || plastic_model == PlasticModel::DPVisc || plastic_model == PlasticModel::MCCVisc || plastic_model == PlasticModel::DPMui || plastic_model == PlasticModel::MCCMui){ Eigen::JacobiSVD svd(Fe_trial, Eigen::ComputeFullU | Eigen::ComputeFullV); // TV hencky = svd.singularValues().array().log(); @@ -26,7 +26,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F hencky_deviatoric /= hencky_deviatoric_norm; // normalize the deviatoric vector so it gives a unit vector specifying the deviatoric direction - if (plastic_model == VM){ + if (plastic_model == PlasticModel::VM){ T q_yield = q_max; @@ -47,7 +47,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F } // end VonMises - else if (plastic_model == DP){ + else if (plastic_model == PlasticModel::DP){ // trial stresses T p_trial = -K * hencky_trace; @@ -83,7 +83,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F } // end DP - else if (plastic_model == DPSoft){ + else if (plastic_model == PlasticModel::DPSoft){ // trial stresses T p_trial = -K * hencky_trace; @@ -160,7 +160,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F } // end DPSoft - else if (plastic_model == VMVisc){ + else if (plastic_model == PlasticModel::VMVisc){ // trial T q_trial = e_mu_prefac * hencky_deviatoric_norm; @@ -236,7 +236,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F } // end VMVisc - else if (plastic_model == DPVisc){ + else if (plastic_model == PlasticModel::DPVisc){ // trial stresses T p_trial = -K * hencky_trace; @@ -338,7 +338,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F } // end DPVisc - else if (plastic_model == DPMui){ + else if (plastic_model == PlasticModel::DPMui){ // trial stresses T p_trial = -K * hencky_trace; @@ -401,7 +401,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F } // end DPMui - else if (plastic_model == MCCMui) { + else if (plastic_model == PlasticModel::MCCMui) { // the trial stress states T p_stress = -K * hencky_trace; @@ -413,22 +413,22 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F ////////////////////////////////////////////////////////////////////// bool perform_rma; T p_c; - if (hardening_law == ExpoExpl){ // Exponential Explicit Hardening + if (hardening_law == HardeningLaw::ExpoExpl){ // Exponential Explicit Hardening p_c = std::max(stress_tolerance, p0 * std::exp(-xi*particles.eps_pl_vol[p])); - perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p_c, beta, mu, K, rma_prefac); - // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p_c, beta, mu, K, rma_prefac); + perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, mu_1, p_c, beta, mu, K, f_mu_prefac); + // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, mu_1, p_c, beta, mu, K); } - else if (hardening_law == SinhExpl){ // Sinh Explicit Hardening + else if (hardening_law == HardeningLaw::SinhExpl){ // Sinh Explicit Hardening p_c = std::max(stress_tolerance, K * std::sinh(-xi*particles.eps_pl_vol[p] + std::asinh(p0/K))); - perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p_c, beta, mu, K, rma_prefac); - // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p_c, beta, mu, K, rma_prefac); + perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, mu_1, p_c, beta, mu, K, f_mu_prefac); + // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, mu_1, p_c, beta, mu, K); } - else if (hardening_law == ExpoImpl){ // Exponential Implicit Hardening - perform_rma = MCCRMAImplicitExponentialOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, rma_prefac, particles.eps_pl_vol[p]); - // perform_rma = MCCRMAImplicitExponential(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, rma_prefac, particles.eps_pl_vol[p]); + else if (hardening_law == HardeningLaw::ExpoImpl){ // Exponential Implicit Hardening + perform_rma = MCCRMAImplicitExponentialOnevar(p_stress, q_stress, exit, mu_1, p0, beta, mu, K, xi, particles.eps_pl_vol[p]); + // perform_rma = MCCRMAImplicitExponential(p_stress, q_stress, exit, mu_1, p0, beta, mu, K, xi, f_mu_prefac, particles.eps_pl_vol[p]); } - else if (hardening_law == SinhImpl) { - perform_rma = MCCRMAImplicitSinhOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, rma_prefac, particles.eps_pl_vol[p]); + else if (hardening_law == HardeningLaw::SinhImpl) { + perform_rma = MCCRMAImplicitSinhOnevar(p_stress, q_stress, exit, mu_1, p0, beta, mu, K, xi, particles.eps_pl_vol[p]); // perform_rma = } else{ @@ -449,10 +449,10 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F particles.eps_pl_vol[p] += eps_pl_vol_inst; if (q_trial > (q_stress + stress_tolerance)) { - if (hardening_law == ExpoImpl){ + if (hardening_law == HardeningLaw::ExpoImpl){ p_c = std::max( p0 * std::exp(-xi*particles.eps_pl_vol[p]) , p_stress + stress_tolerance ); } - else if (hardening_law == SinhImpl){ + else if (hardening_law == HardeningLaw::SinhImpl){ p_c = std::max( K * std::sinh(-xi*particles.eps_pl_vol[p] + std::asinh(p0/K)) , p_stress + stress_tolerance ); } T p_special = p_stress + beta * p_c + stress_tolerance; // always larger than 0 @@ -483,7 +483,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F } // end MCCMui - else if (plastic_model == MCC || plastic_model == MCCVisc){ + else if (plastic_model == PlasticModel::MCC || plastic_model == PlasticModel::MCCVisc){ // the trial stress states T p_stress = -K * hencky_trace; @@ -493,33 +493,33 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F T p_trial = p_stress; T q_trial = q_stress; - if (plastic_model == MCCVisc){ + if (plastic_model == PlasticModel::MCCVisc){ if (use_mibf) particles.muI[p] = M; } bool perform_rma; T p_c; - if (hardening_law == NoHard){ // Exponential Explicit Hardening - perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p0, beta, mu, K, rma_prefac); - // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K, rma_prefac); + if (hardening_law == HardeningLaw::NoHard){ // Exponential Explicit Hardening + perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p0, beta, mu, K, f_mu_prefac); + // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K); } - else if (hardening_law == ExpoExpl){ // Exponential Explicit Hardening + else if (hardening_law == HardeningLaw::ExpoExpl){ // Exponential Explicit Hardening p_c = std::max(stress_tolerance, p0*std::exp(-xi*particles.eps_pl_vol[p])); - perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p_c, beta, mu, K, rma_prefac); - // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p_c, beta, mu, K, rma_prefac); + perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p_c, beta, mu, K, f_mu_prefac); + // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p_c, beta, mu, K); } - else if (hardening_law == SinhExpl){ // Sinh Explicit Hardening + else if (hardening_law == HardeningLaw::SinhExpl){ // Sinh Explicit Hardening p_c = std::max(stress_tolerance, K*std::sinh(-xi*particles.eps_pl_vol[p] + std::asinh(p0/K))); - perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p_c, beta, mu, K, rma_prefac); - // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p_c, beta, mu, K, rma_prefac); + perform_rma = MCCRMAExplicit(p_stress, q_stress, exit, M, p_c, beta, mu, K, f_mu_prefac); + // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, M, p_c, beta, mu, K); } - else if (hardening_law == ExpoImpl){ // Exponential Implicit Hardening - perform_rma = MCCRMAImplicitExponentialOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, rma_prefac, particles.eps_pl_vol[p]); - // perform_rma = MCCRMAImplicitExponential(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, rma_prefac, particles.eps_pl_vol[p]); + else if (hardening_law == HardeningLaw::ExpoImpl){ // Exponential Implicit Hardening + perform_rma = MCCRMAImplicitExponentialOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, particles.eps_pl_vol[p]); + // perform_rma = MCCRMAImplicitExponential(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, f_mu_prefac, particles.eps_pl_vol[p]); } - else if (hardening_law == SinhImpl) { - perform_rma = MCCRMAImplicitSinhOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, rma_prefac, particles.eps_pl_vol[p]); + else if (hardening_law == HardeningLaw::SinhImpl) { + perform_rma = MCCRMAImplicitSinhOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, particles.eps_pl_vol[p]); } else{ debug("You specified an invalid HARDENING LAW!"); @@ -536,7 +536,7 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F T delta_gamma; //////////////////////////////////////////////////////////////// - if (plastic_model == MCCVisc){ + if (plastic_model == PlasticModel::MCCVisc){ T q_yield = q_stress; @@ -591,10 +591,10 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F q_stress = std::max(q_yield, q_trial - f_mu_prefac * delta_gamma); // delta_gamma = dt * gamma_dot_S if (use_mibf){ - if (hardening_law == ExpoImpl){ + if (hardening_law == HardeningLaw::ExpoImpl){ p_c = std::max(T(0), p0 * std::exp(-xi*particles.eps_pl_vol[p])); } - else if (hardening_law == SinhImpl){ + else if (hardening_law == HardeningLaw::SinhImpl){ p_c = std::max(T(0), K * std::sinh(-xi*particles.eps_pl_vol[p] + std::asinh(p0/K))); } diff --git a/src/simulation/position_update.cpp b/src/simulation/position_update.cpp index caa6aef..b921b80 100644 --- a/src/simulation/position_update.cpp +++ b/src/simulation/position_update.cpp @@ -5,7 +5,7 @@ void Simulation::positionUpdate(){ - #pragma omp parallel for num_threads(n_threads) + #pragma omp parallel for schedule(static) num_threads(n_threads) for(int p=0; pmove(dt, frame_dt, time); } } ///////////////////////////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////// EXTRA HELPER FUNCTIONS ////////////////////////////////////////////////// +////////////////////////////////////// EXTRA HELPER FUNCTIONS /////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////////////////////////// T Simulation::calculateBulkModulus(){ diff --git a/src/simulation/simulation.hpp b/src/simulation/simulation.hpp index 432cca1..9a3acac 100644 --- a/src/simulation/simulation.hpp +++ b/src/simulation/simulation.hpp @@ -11,6 +11,7 @@ #include #include #include +#include #include "../tools.hpp" #include "../data_structures.hpp" @@ -19,7 +20,6 @@ #include "../objects/object_general.hpp" #include "../objects/object_plate.hpp" - class Simulation{ public: Simulation(); @@ -27,9 +27,14 @@ class Simulation{ int exit = 0; +#ifdef THREEDIM + const unsigned int dim = 3; +#else + const unsigned int dim = 2; +#endif + unsigned int n_threads = 1; unsigned int end_frame = 1; - T fps = 1; bool is_initialized = false; bool save_sim = true; @@ -41,47 +46,44 @@ class Simulation{ bool use_mibf = false; bool use_musl = false; + TV grid_reference_point = 2e10 * TV::Ones(); + TV gravity = TV::Zero(); + + T fps = 1; T cfl = 0.5; T cfl_elastic = 0.5; T flip_ratio = -0.95; - T rho = 1000; - - TV gravity = TV::Zero(); T gravity_time = 0; // bool no_liftoff = true; - T Lx = 1; T Ly = 1; #ifdef THREEDIM T Lz = 1; #endif - TV grid_reference_point = 2e10 * TV::Ones(); - // Particle data Particles particles; unsigned int Np; + unsigned int num_add_pbc_particles; T particle_mass; T particle_volume; // initial particle volume T dx; // Elastoplasticity - ElasticModel elastic_model = Hencky; - PlasticModel plastic_model = NoPlasticity; - HardeningLaw hardening_law = ExpoImpl; - + ElasticModel elastic_model = ElasticModel::Hencky; + PlasticModel plastic_model = PlasticModel::NoPlasticity; + HardeningLaw hardening_law = HardeningLaw::ExpoImpl; + bool use_pradhana = true; T E = 1e6; // Young's modulus (3D) T nu = 0.3; // Poisson's ratio (3D) - - bool use_pradhana = true; - bool use_mises_q = false; + T stress_tolerance = 1e-5; + T xi = 0; // Von Mises: T q_max = 100; T q_min = 100; T p_min = -1e20; - T xi = 0; // Drucker Prager T M = 1; @@ -102,11 +104,12 @@ class Simulation{ T mu_1 = std::tan(20.9*M_PI/180.0); T mu_2 = std::tan(32.8*M_PI/180.0);; - T stress_tolerance = 1e-5; + // Prefactor for q in plasticity models + T q_prefac = 1.0 / std::sqrt(2.0); // q = factor * ||dev(tau)|| // Objects - std::vector plates; - std::vector objects; + std::vector> plates; + std::vector> objects; // Functions void initialize(bool save = true, std::string dir = "output/", std::string name = "dummy"); @@ -117,59 +120,39 @@ class Simulation{ void computeAvgData(TM& volavg_cauchy, TM& volavg_kirchh, T& Javg); void saveParticleData(std::string extra = ""); void saveGridData(std::string extra = ""); - void createDirectory(); - void advanceStep(); void updateDt(); - void resizeGrid(); void remeshFixed(unsigned int extra_nodes); void remeshFixedInit(unsigned int sfx, unsigned int sfy, unsigned int sfz); void remeshFixedCont(); - void P2G(); void explicitEulerUpdate(); void G2P(); void deformationUpdate(); - void MUSL(); - void positionUpdate(); void PBCAddParticles1D(); void PBCAddParticles(unsigned int safety_factor); void PBCDelParticles(); - unsigned int num_add_pbc_particles; - - TM NeoHookeanPiola(TM & Fe); - TM HenckyPiola(TM & Fe); void plasticity(unsigned int p, unsigned int & plastic_count, TM & Fe_trial); - void moveObjects(); - void boundaryCollision(int index, TV Xi, TV& vi); void overwriteGridVelocity(TV Xi, TV& vi); - - T calculateBulkModulus(); void checkMomentumConservation(); void checkMassConservation(); - void addExternalParticleGravity(); std::pair createExternalGridGravity(); - - - #ifdef THREEDIM - const unsigned int dim = 3; - #else - const unsigned int dim = 2; - #endif + T calculateBulkModulus(); + TM NeoHookeanPiola(TM & Fe); + TM HenckyPiola(TM & Fe); private: unsigned int current_time_step = 0; unsigned int frame = 0; T time = 0; - T runtime_p2g = 0; T runtime_g2p = 0; T runtime_euler = 0; @@ -179,26 +162,22 @@ class Simulation{ std::string sim_name; std::string directory; + TV gravity_final; + T final_time; T frame_dt; T dt; T dt_max; T wave_speed; - - TV gravity_final; - T mu; // shear modulus T lambda; // first Lame parameter T K; // bulk modulus - T fac_Q; // for mu(I) rheology // Prefactors for plasticity models - T q_prefac; // q = factor * ||dev(tau)|| T d_prefac; // gamma = factor * ||dev(eps)|| T e_mu_prefac; // q = factor * ||dev(eps)|| T f_mu_prefac; // q^tr - q = factor * dt * gamma_dot - T rma_prefac; // Precomputations T one_over_dx; @@ -214,13 +193,13 @@ class Simulation{ unsigned int grid_nodes; #ifdef THREEDIM - inline unsigned int ind(unsigned int i, unsigned int j, unsigned int k){ - return (i*Ny + j) * Nz + k; // 3D - } + inline unsigned int ind(unsigned int i, unsigned int j, unsigned int k){ + return (i*Ny + j) * Nz + k; // 3D + } #else - inline unsigned int ind(unsigned int i, unsigned int j){ - return (i*Ny + j); // 3D - } + inline unsigned int ind(unsigned int i, unsigned int j){ + return (i*Ny + j); // 3D + } #endif T min_x_init; @@ -242,11 +221,7 @@ class Simulation{ T low_z_init; T high_z_init; #endif - - -}; // END Simulation class - - +}; // end Simulation class inline TM Simulation::NeoHookeanPiola(TM & Fe){ return mu * (Fe - Fe.transpose().inverse()) + lambda * std::log(Fe.determinant()) * Fe.transpose().inverse(); @@ -261,6 +236,4 @@ inline TM Simulation::HenckyPiola(TM & Fe){ return dPsidF; } // end HenckyPiola - - #endif // SIMULATION_HPP diff --git a/src/simulation/update_dt.cpp b/src/simulation/update_dt.cpp index 69b2484..30f7e44 100644 --- a/src/simulation/update_dt.cpp +++ b/src/simulation/update_dt.cpp @@ -14,7 +14,6 @@ void Simulation::updateDt(){ if (max_speed >= wave_speed){ debug(" Detected particle speed ", max_speed, " larger than elastic wave speed ", wave_speed); // exit = 1; - return; } #ifdef WARNINGS @@ -42,6 +41,8 @@ void Simulation::updateDt(){ #endif + // Here one may hard-code a special gravity evolution with time + // Default: linear gravity increase until "gravity_time" if (gravity_special){ if (time < gravity_time){ diff --git a/src/tools.hpp b/src/tools.hpp index 86d04e9..14be3a0 100644 --- a/src/tools.hpp +++ b/src/tools.hpp @@ -50,11 +50,11 @@ typedef double T; // float or double #endif //////////////////////// -enum PlateType { top, bottom, left, right, front, back }; -enum ElasticModel { Hencky, NeoHookean }; -enum PlasticModel { NoPlasticity, VM, DP, DPSoft, MCC, VMVisc, DPVisc, MCCVisc, DPMui, MCCMui}; -enum HardeningLaw { NoHard, ExpoExpl, ExpoImpl, SinhExpl, SinhImpl }; -enum BoundaryCondition { NOSLIP, SLIPSTICK, SLIPFREE }; +enum class PlateType { top, bottom, left, right, front, back }; +enum class ElasticModel { Hencky, NeoHookean }; +enum class PlasticModel { NoPlasticity, VM, DP, DPSoft, MCC, VMVisc, DPVisc, MCCVisc, DPMui, MCCMui}; +enum class HardeningLaw { NoHard, ExpoExpl, ExpoImpl, SinhExpl, SinhImpl }; +enum class BC { NoSlip, SlipStick, SlipFree }; ///////////////////// TOOLS //////////////////////// diff --git a/tests/tests.cpp b/tests/tests.cpp index deb0400..07c3bcf 100644 --- a/tests/tests.cpp +++ b/tests/tests.cpp @@ -10,10 +10,11 @@ #include "../src/objects/object_ground.hpp" -TEST(BoundaryTest, AnalyticSLIPFREE) { +TEST(BoundaryTest, AnalyticSlipFree) { Simulation sim; sim.initialize(false); + sim.reduce_verbose = true; T half_rounds = 1; int num_frames_in_half_round = 100; @@ -24,8 +25,8 @@ TEST(BoundaryTest, AnalyticSLIPFREE) { sim.cfl = 0.5; sim.flip_ratio = -0.95; - sim.elastic_model = Hencky; - sim.plastic_model = NoPlasticity; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::NoPlasticity; sim.E = 1e6; sim.nu = 0.3; sim.rho = 1000; @@ -41,21 +42,23 @@ TEST(BoundaryTest, AnalyticSLIPFREE) { sim.particles.x[0](0) = -1; sim.particles.x[0](1) = 1; - ObjectCurve curve = ObjectCurve(SLIPFREE, 0.0); sim.objects.push_back(&curve); + sim.objects.push_back(std::make_unique(BC::SlipFree, 0.0)); sim.simulate(); T v_sim = sim.particles.v[0](0); T v_true = std::sqrt(2*9.81); T diff = std::abs(v_sim-v_true)/v_true; - ASSERT_NEAR(diff, 0.0, 1e-4); + debug("diff: ", diff); + ASSERT_NEAR(diff, 0.0, 1e-3); } -TEST(BoundaryTest, AnalyticSLIPSTICK) { +TEST(BoundaryTest, AnalyticSlipStick) { Simulation sim; sim.initialize(false); + sim.reduce_verbose = true; T half_rounds = 1; int num_frames_in_half_round = 100; @@ -66,8 +69,8 @@ TEST(BoundaryTest, AnalyticSLIPSTICK) { sim.cfl = 0.5; sim.flip_ratio = -0.95; - sim.elastic_model = Hencky; - sim.plastic_model = NoPlasticity; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::NoPlasticity; sim.E = 1e6; sim.nu = 0.3; sim.rho = 1000; @@ -81,22 +84,24 @@ TEST(BoundaryTest, AnalyticSLIPSTICK) { sim.particles = Particles(sim.Np); sim.particles.x[0](0) = -1; - sim.particles.x[0](1) = 0.9999; + sim.particles.x[0](1) = 1; - ObjectCurve curve = ObjectCurve(SLIPSTICK, 0.0); sim.objects.push_back(&curve); + sim.objects.push_back(std::make_unique(BC::SlipStick, 0.0)); sim.simulate(); T v_sim = sim.particles.v[0](0); T v_true = std::sqrt(2*9.81); T diff = std::abs(v_sim-v_true)/v_true; + debug("diff: ", diff); ASSERT_NEAR(diff, 0.0, 1e-3); } -TEST(CoulombFrictionTest, PlateSLIPFREE) { +TEST(CoulombFrictionTest, PlateSlipFree) { Simulation sim; sim.initialize(false); + sim.reduce_verbose = true; sim.end_frame = 5; sim.fps = 1; @@ -110,8 +115,8 @@ TEST(CoulombFrictionTest, PlateSLIPFREE) { sim.gravity[0] = +9.81 * std::sin(theta); sim.gravity[1] = -9.81 * std::cos(theta); - sim.elastic_model = Hencky; - sim.plastic_model = NoPlasticity; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::NoPlasticity; sim.E = 1e6; sim.nu = 0.3; sim.rho = 1000; @@ -128,8 +133,7 @@ TEST(CoulombFrictionTest, PlateSLIPFREE) { T friction = std::tan(15.0 * M_PI / 180.0); - ObjectPlate ground = ObjectPlate(0, bottom, SLIPFREE, friction); - sim.plates.push_back(ground); + sim.plates.push_back(std::make_unique(0, PlateType::bottom, BC::SlipFree, friction)); sim.simulate(); @@ -141,13 +145,15 @@ TEST(CoulombFrictionTest, PlateSLIPFREE) { T final_time = (T)sim.end_frame / sim.fps; T mean_true_x = 0.5*9.81*final_time*final_time*(std::sin(theta) - friction*std::cos(theta)); T diff = std::abs(mean_sim_x-mean_true_x) / mean_true_x; + debug("diff: ", diff); ASSERT_NEAR(diff, 0.0, 1e-3); } -TEST(CoulombFrictionTest, PlateSLIPSTICK) { +TEST(CoulombFrictionTest, PlateSlipStick) { Simulation sim; sim.initialize(false); + sim.reduce_verbose = true; sim.end_frame = 5; sim.fps = 1; @@ -161,8 +167,8 @@ TEST(CoulombFrictionTest, PlateSLIPSTICK) { sim.gravity[0] = +9.81 * std::sin(theta); sim.gravity[1] = -9.81 * std::cos(theta); - sim.elastic_model = Hencky; - sim.plastic_model = NoPlasticity; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::NoPlasticity; sim.E = 1e6; sim.nu = 0.3; sim.rho = 1000; @@ -179,8 +185,7 @@ TEST(CoulombFrictionTest, PlateSLIPSTICK) { T friction = std::tan(15.0 * M_PI / 180.0); - ObjectPlate ground = ObjectPlate(0, bottom, SLIPSTICK, friction); - sim.plates.push_back(ground); + sim.plates.push_back(std::make_unique(0, PlateType::bottom, BC::SlipStick, friction)); sim.simulate(); @@ -192,15 +197,15 @@ TEST(CoulombFrictionTest, PlateSLIPSTICK) { T final_time = (T)sim.end_frame / sim.fps; T mean_true_x = 0.5*9.81*final_time*final_time*(std::sin(theta) - friction*std::cos(theta)); T diff = std::abs(mean_sim_x-mean_true_x) / mean_true_x; + debug("diff: ", diff); ASSERT_NEAR(diff, 0.0, 1e-3); } - - -TEST(CoulombFrictionTest, GeneralSLIPFREE) { +TEST(CoulombFrictionTest, GeneralSlipFree) { Simulation sim; sim.initialize(false); + sim.reduce_verbose = true; sim.end_frame = 5; sim.fps = 1; @@ -214,8 +219,8 @@ TEST(CoulombFrictionTest, GeneralSLIPFREE) { sim.gravity[0] = +9.81 * std::sin(theta); sim.gravity[1] = -9.81 * std::cos(theta); - sim.elastic_model = Hencky; - sim.plastic_model = NoPlasticity; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::NoPlasticity; sim.E = 1e6; sim.nu = 0.3; sim.rho = 1000; @@ -232,8 +237,7 @@ TEST(CoulombFrictionTest, GeneralSLIPFREE) { T friction = std::tan(15.0 * M_PI / 180.0); - ObjectGround ground = ObjectGround(SLIPFREE, friction); - sim.objects.push_back(&ground); + sim.objects.push_back(std::make_unique(BC::SlipFree, friction)); sim.simulate(); @@ -245,13 +249,15 @@ TEST(CoulombFrictionTest, GeneralSLIPFREE) { T final_time = (T)sim.end_frame / sim.fps; T mean_true_x = 0.5*9.81*final_time*final_time*(std::sin(theta) - friction*std::cos(theta)); T diff = std::abs(mean_sim_x-mean_true_x) / mean_true_x; + debug("diff: ", diff); ASSERT_NEAR(diff, 0.0, 1e-3); } -TEST(CoulombFrictionTest, GeneralSLIPSTICK) { +TEST(CoulombFrictionTest, GeneralSlipStick) { Simulation sim; sim.initialize(false); + sim.reduce_verbose = true; sim.end_frame = 5; sim.fps = 1; @@ -265,8 +271,8 @@ TEST(CoulombFrictionTest, GeneralSLIPSTICK) { sim.gravity[0] = +9.81 * std::sin(theta); sim.gravity[1] = -9.81 * std::cos(theta); - sim.elastic_model = Hencky; - sim.plastic_model = NoPlasticity; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::NoPlasticity; sim.E = 1e6; sim.nu = 0.3; sim.rho = 1000; @@ -283,8 +289,7 @@ TEST(CoulombFrictionTest, GeneralSLIPSTICK) { T friction = std::tan(15.0 * M_PI / 180.0); - ObjectGround ground = ObjectGround(SLIPSTICK, friction); - sim.objects.push_back(&ground); + sim.objects.push_back(std::make_unique(BC::SlipStick, friction)); sim.simulate(); @@ -296,19 +301,15 @@ TEST(CoulombFrictionTest, GeneralSLIPSTICK) { T final_time = (T)sim.end_frame / sim.fps; T mean_true_x = 0.5*9.81*final_time*final_time*(std::sin(theta) - friction*std::cos(theta)); T diff = std::abs(mean_sim_x-mean_true_x) / mean_true_x; + debug("diff: ", diff); ASSERT_NEAR(diff, 0.0, 1e-3); } - - - TEST(BoundaryTest, MIBF) { T friction = std::tan(30.0 * M_PI / 180.0); T theta = 32 * M_PI / 180; - - ObjectPlate ground = ObjectPlate(0, bottom, SLIPFREE, friction); - + Simulation sim_one; sim_one.initialize(false); @@ -323,7 +324,7 @@ TEST(BoundaryTest, MIBF) { sim_one.gravity[0] = +9.81 * std::sin(theta); sim_one.gravity[1] = -9.81 * std::cos(theta); - sim_one.elastic_model = Hencky; + sim_one.elastic_model = ElasticModel::Hencky; sim_one.E = 1e5; sim_one.nu = 0.3; sim_one.rho = 1000; @@ -340,13 +341,12 @@ TEST(BoundaryTest, MIBF) { } sim_one.grid_reference_point = TV::Zero(); - sim_one.plates.push_back(ground); + sim_one.plates.push_back(std::make_unique(0, PlateType::bottom, BC::SlipFree, friction)); - sim_one.plastic_model = DPVisc; + sim_one.plastic_model = PlasticModel::DPVisc; sim_one.use_pradhana = false; - sim_one.use_mises_q = false; - sim_one.use_mibf = false; // NB + sim_one.use_mibf = false; // NB sim_one.M = friction; sim_one.q_cohesion = 0; @@ -375,7 +375,7 @@ TEST(BoundaryTest, MIBF) { sim_two.gravity[0] = +9.81 * std::sin(theta); sim_two.gravity[1] = -9.81 * std::cos(theta); - sim_two.elastic_model = Hencky; + sim_two.elastic_model = ElasticModel::Hencky; sim_two.E = 1e5; sim_two.nu = 0.3; sim_two.rho = 1000; @@ -392,13 +392,12 @@ TEST(BoundaryTest, MIBF) { } sim_two.grid_reference_point = TV::Zero(); - sim_two.plates.push_back(ground); + sim_two.plates.push_back(std::make_unique(0, PlateType::bottom, BC::SlipFree, friction)); - sim_two.plastic_model = DPVisc; + sim_two.plastic_model = PlasticModel::DPVisc; sim_two.use_pradhana = false; - sim_two.use_mises_q = false; - sim_two.use_mibf = true; // NB + sim_two.use_mibf = true; // NB sim_two.M = friction; sim_two.q_cohesion = 0; @@ -415,6 +414,8 @@ TEST(BoundaryTest, MIBF) { ASSERT_NEAR(diff, 0.0, 1e-13); + + Simulation sim_three; sim_three.initialize(false); @@ -429,7 +430,7 @@ TEST(BoundaryTest, MIBF) { sim_three.gravity[0] = +9.81 * std::sin(theta); sim_three.gravity[1] = -9.81 * std::cos(theta); - sim_three.elastic_model = Hencky; + sim_three.elastic_model = ElasticModel::Hencky; sim_three.E = 1e5; sim_three.nu = 0.3; sim_three.rho = 1000; @@ -446,15 +447,16 @@ TEST(BoundaryTest, MIBF) { } sim_three.grid_reference_point = TV::Zero(); - sim_three.plates.push_back(ground); + sim_three.plates.push_back(std::make_unique(0, PlateType::bottom, BC::SlipFree, friction)); - sim_three.plastic_model = DPVisc; + sim_three.plastic_model = PlasticModel::DPVisc; sim_three.use_pradhana = false; - sim_three.use_mises_q = true; // NB - sim_three.use_mibf = false; // NB + sim_three.use_mibf = false; // NB + sim_three.q_prefac = std::sqrt(3.0/2.0); // NB - sim_three.M = std::sqrt(3.0)*friction; // NB + + sim_three.M = std::sqrt(3.0)*friction; // NB sim_three.q_cohesion = 0; sim_three.perzyna_exp = 1; sim_three.perzyna_visc = 0; @@ -470,8 +472,6 @@ TEST(BoundaryTest, MIBF) { } - - TEST(ElasticityTest, BulkModulus) { Simulation sim; @@ -502,15 +502,15 @@ TEST(ElasticityTest, BulkModulus) { T load_factor = 1; #ifdef THREEDIM - ObjectPlate ground = ObjectPlate(0-0.5*sim.dx, bottom, NOSLIP, 0, -1e15, 1e15, 0, vel, 0, vmin_factor, load_factor); sim.plates.push_back(ground); - ObjectPlate compre = ObjectPlate(sim.Ly+0.5*sim.dx, top, NOSLIP, 0, -1e15, 1e15, 0, -vel, 0, vmin_factor, load_factor); sim.plates.push_back(compre); + sim.plates.push_back(std::make_unique(0-0.5*sim.dx, PlateType::bottom, BC::NoSlip, 0, -1e15, 1e15, 0, vel, 0, vmin_factor, load_factor)); + sim.plates.push_back(std::make_unique(sim.Ly+0.5*sim.dx, PlateType::top, BC::NoSlip, 0, -1e15, 1e15, 0, -vel, 0, vmin_factor, load_factor)); #else - ObjectPlate ground = ObjectPlate(0-0.5*sim.dx, bottom, NOSLIP, 0, -1e15, 1e15, 0, vel, vmin_factor, load_factor); sim.plates.push_back(ground); - ObjectPlate compre = ObjectPlate(sim.Ly+0.5*sim.dx, top, NOSLIP, 0, -1e15, 1e15, 0, -vel, vmin_factor, load_factor); sim.plates.push_back(compre); + sim.plates.push_back(std::make_unique(0-0.5*sim.dx, PlateType::bottom, BC::NoSlip, 0, -1e15, 1e15, 0, vel, vmin_factor, load_factor)); + sim.plates.push_back(std::make_unique(sim.Ly+0.5*sim.dx, PlateType::top, BC::NoSlip, 0, -1e15, 1e15, 0, -vel, vmin_factor, load_factor)); #endif - sim.elastic_model = Hencky; - sim.plastic_model = NoPlasticity; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::NoPlasticity; sim.simulate(); @@ -532,12 +532,13 @@ TEST(ElasticityTest, BulkModulus) { T rel_diff = std::abs(measured_K - true_K) / true_K; + debug("true_K: ", true_K); + debug("measured_K: ", measured_K); + debug("rel_diff: ", rel_diff); + ASSERT_NEAR(rel_diff, 0.0, 0.03); } - - - TEST(EnergyTest, Rotation) { Simulation sim; @@ -574,8 +575,8 @@ TEST(EnergyTest, Rotation) { total_energy_init += 0.5*(vx*vx + vy*vy); // per unit mass } - sim.elastic_model = Hencky; - sim.plastic_model = NoPlasticity; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::NoPlasticity; sim.simulate(); @@ -593,14 +594,13 @@ TEST(EnergyTest, Rotation) { T rel_diff = (total_energy_init - total_energy_last) / total_energy_init; + debug("rel_diff: ", rel_diff); + EXPECT_TRUE((rel_diff >= 0) && (rel_diff <= 1e-3)); // Can not have energy increase! // Energy decrease within a relative tolerance } - - - TEST(CollapseTest, DruckerPragerOne) { Simulation sim; @@ -637,10 +637,9 @@ TEST(CollapseTest, DruckerPragerOne) { } sim.grid_reference_point = TV::Zero(); - sim.elastic_model = Hencky; - sim.plastic_model = DPSoft; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::DPSoft; - sim.use_mises_q = false; sim.use_pradhana = true; sim.xi = 0; @@ -649,13 +648,13 @@ TEST(CollapseTest, DruckerPragerOne) { sim.M = std::tan(30.0 * M_PI / 180.0); #ifdef THREEDIM - ObjectPlate sideback = ObjectPlate(0, back, SLIPFREE); sim.plates.push_back(sideback); - ObjectPlate sidefront = ObjectPlate(sim.Lz, front, SLIPFREE); sim.plates.push_back(sidefront); - ObjectPlate sideleft = ObjectPlate(0, left, SLIPFREE); sim.plates.push_back(sideleft); - ObjectPlate ground = ObjectPlate(0, bottom, NOSLIP); sim.plates.push_back(ground); + sim.plates.push_back(std::make_unique(0, PlateType::back, BC::SlipFree)); + sim.plates.push_back(std::make_unique(sim.Lz, PlateType::front, BC::SlipFree)); + sim.plates.push_back(std::make_unique(0, PlateType::left, BC::SlipFree)); + sim.plates.push_back(std::make_unique(0, PlateType::bottom, BC::NoSlip)); #else - ObjectPlate sideleft = ObjectPlate(0, left, SLIPFREE); sim.plates.push_back(sideleft); - ObjectPlate ground = ObjectPlate(0, bottom, NOSLIP); sim.plates.push_back(ground); + sim.plates.push_back(std::make_unique(0, PlateType::left, BC::SlipFree)); + sim.plates.push_back(std::make_unique(0, PlateType::bottom, BC::NoSlip)); #endif sim.simulate(); @@ -663,6 +662,7 @@ TEST(CollapseTest, DruckerPragerOne) { auto max_x_it = std::max_element( sim.particles.x.begin(), sim.particles.x.end(), [](const TV &x1, const TV &x2){return x1(0) < x2(0);} ); T max_x = (*max_x_it)(0); T diff = std::abs(max_x - 0.56); + debug("diff: ", diff); ASSERT_NEAR(diff, 0.0, 0.011); } @@ -702,10 +702,9 @@ TEST(CollapseTest, DruckerPragerTwo) { } sim.grid_reference_point = TV::Zero(); - sim.elastic_model = Hencky; - sim.plastic_model = DPVisc; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::DPVisc; - sim.use_mises_q = false; sim.use_pradhana = true; sim.perzyna_visc = 0; @@ -715,13 +714,13 @@ TEST(CollapseTest, DruckerPragerTwo) { sim.M = std::tan(30.0 * M_PI / 180.0); #ifdef THREEDIM - ObjectPlate sideback = ObjectPlate(0, back, SLIPFREE); sim.plates.push_back(sideback); - ObjectPlate sidefront = ObjectPlate(sim.Lz, front, SLIPFREE); sim.plates.push_back(sidefront); - ObjectPlate sideleft = ObjectPlate(0, left, SLIPFREE); sim.plates.push_back(sideleft); - ObjectPlate ground = ObjectPlate(0, bottom, NOSLIP); sim.plates.push_back(ground); + sim.plates.push_back(std::make_unique(0, PlateType::back, BC::SlipFree)); + sim.plates.push_back(std::make_unique(sim.Lz, PlateType::front, BC::SlipFree)); + sim.plates.push_back(std::make_unique(0, PlateType::left, BC::SlipFree)); + sim.plates.push_back(std::make_unique(0, PlateType::bottom, BC::NoSlip)); #else - ObjectPlate sideleft = ObjectPlate(0, left, SLIPFREE); sim.plates.push_back(sideleft); - ObjectPlate ground = ObjectPlate(0, bottom, NOSLIP); sim.plates.push_back(ground); + sim.plates.push_back(std::make_unique(0, PlateType::left, BC::SlipFree)); + sim.plates.push_back(std::make_unique(0, PlateType::bottom, BC::NoSlip)); #endif sim.simulate(); @@ -729,5 +728,6 @@ TEST(CollapseTest, DruckerPragerTwo) { auto max_x_it = std::max_element( sim.particles.x.begin(), sim.particles.x.end(), [](const TV &x1, const TV &x2){return x1(0) < x2(0);} ); T max_x = (*max_x_it)(0); T diff = std::abs(max_x - 0.56); + debug("diff: ", diff); ASSERT_NEAR(diff, 0.0, 0.011); } diff --git a/tests/tests_vdb.cpp b/tests/tests_vdb.cpp index 3859a35..6287d2c 100644 --- a/tests/tests_vdb.cpp +++ b/tests/tests_vdb.cpp @@ -14,7 +14,6 @@ TEST(FileIOTest, LoadFile) { ASSERT_TRUE(file.is_open()); } - TEST(BoundaryTest, CalcNormal) { openvdb::initialize(); @@ -36,10 +35,8 @@ TEST(BoundaryTest, CalcNormal) { EXPECT_TRUE(obj1.inside(pos) == obj2.inside(pos)); ASSERT_NEAR(diff, 0.0, 1e-4); - } - TEST(BoundaryTest, VDB) { openvdb::initialize(); @@ -56,8 +53,8 @@ TEST(BoundaryTest, VDB) { sim.cfl = 0.5; sim.flip_ratio = -0.95; - sim.elastic_model = Hencky; - sim.plastic_model = NoPlasticity; + sim.elastic_model = ElasticModel::Hencky; + sim.plastic_model = PlasticModel::NoPlasticity; sim.E = 1e6; sim.nu = 0.3; sim.rho = 1000; @@ -74,8 +71,7 @@ TEST(BoundaryTest, VDB) { sim.particles.x[0](1) = 1; std::string file_path = std::string(INCLUDE_DIR) + "/curve.vdb"; - ObjectVdb curve = ObjectVdb(file_path, SLIPFREE, 0.0); sim.objects.push_back(&curve); - + sim.objects.push_back(std::make_unique(file_path, BC::SlipFree, 0.0)); sim.simulate();