Keeping DIMMatter up to date with master

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 <number of cores for compilation>` 
 
-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 <name of test>`
+
+
 ### 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<ObjectPlate>(
+        /*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)    
 

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<ObjectPlate>(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<ObjectBump>(BC::SlipFree, friction));
+    // sim.objects.push_back(std::make_unique<ObjectGate>(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<ObjectVdb>("../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

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)

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<TV> x;
@@ -34,12 +32,10 @@ class Particles{
   std::vector<T> eps_pl_dev;
   std::vector<T> eps_pl_vol;
   std::vector<T> eps_pl_vol_pradhana;
-
   std::vector<T> delta_gamma;
   std::vector<T> viscosity;
   std::vector<T> muI;
 
-  std::vector<TM> tau;
   std::vector<TM> F;
   std::vector<TM> Bmat;
 

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<ObjectPlate>(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<ObjectBump>(BC::SlipFree, friction));
+    // sim.objects.push_back(std::make_unique<ObjectGate>(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<ObjectVdb>("../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

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);