Plasticity improvements

README.md

@@ -208,14 +208,14 @@ 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`        | 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:

examples/collapse.cpp

@@ -75,7 +75,7 @@ int main(){
     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/mpm.cpp

@@ -75,7 +75,7 @@ int main(){
     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/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;

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
  */
 {

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;

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
  */
 {

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
  */
 {

src/sampling/sampling_particles.hpp

@@ -122,7 +122,7 @@
 
     } // end sampleParticles
 
-#endif // DIMENSION
+#endif 
 
 
 #endif  // SAMPLING_PARTICLES_HPP

src/sampling/sampling_particles_vdb.hpp

@@ -14,7 +14,7 @@ template <typename S>
 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);

src/simulation/explicit_euler_update.cpp

@@ -16,7 +16,7 @@ void Simulation::explicitEulerUpdate(){
     {
         std::vector<TV> 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];

src/simulation/g2p.cpp

@@ -20,7 +20,7 @@ void Simulation::G2P(){
         std::vector<TV> particles_flip_local(Np); std::fill( particles_flip_local.begin(), particles_flip_local.end(), TV::Zero() );
         std::vector<TM> 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();

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<Np; p++){
 
         //// Velicity is updated
@@ -31,7 +31,7 @@ void Simulation::MUSL(){
     {
         std::vector<TV> 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

src/simulation/p2g.cpp

@@ -15,7 +15,7 @@ void Simulation::P2G(){
         std::vector<T> grid_mass_local(grid_nodes);
         std::vector<T> 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

src/simulation/plasticity.cpp

@@ -415,20 +415,20 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F
             T p_c;
             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, mu_1, p_c, beta, mu, K, rma_prefac);
-                // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, mu_1, 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 == 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, mu_1, p_c, beta, mu, K, rma_prefac);
-                // perform_rma = MCCRMAExplicitOnevar(p_stress, q_stress, exit, mu_1, 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 == HardeningLaw::ExpoImpl){ // Exponential Implicit Hardening
-                   perform_rma = MCCRMAImplicitExponentialOnevar(p_stress, q_stress, exit, mu_1, p0, beta, mu, K, xi, rma_prefac, particles.eps_pl_vol[p]);
-                // perform_rma =       MCCRMAImplicitExponential(p_stress, q_stress, exit, mu_1, p0, beta, mu, K, xi, rma_prefac, particles.eps_pl_vol[p]);
+                   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 == HardeningLaw::SinhImpl) {
-                   perform_rma = MCCRMAImplicitSinhOnevar(p_stress, q_stress, exit, mu_1, p0, beta, mu, K, xi, rma_prefac, particles.eps_pl_vol[p]);
+                   perform_rma = MCCRMAImplicitSinhOnevar(p_stress, q_stress, exit, mu_1, p0, beta, mu, K, xi, particles.eps_pl_vol[p]);
                 // perform_rma = <no alternative at the moment>
             }
             else{
@@ -501,25 +501,25 @@ void Simulation::plasticity(unsigned int p, unsigned int & plastic_count, TM & F
             bool perform_rma;
             T p_c;
             if (hardening_law == HardeningLaw::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);
+                   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 == 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 == 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 == HardeningLaw::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]);
+                   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 == HardeningLaw::SinhImpl) {
-                perform_rma = MCCRMAImplicitSinhOnevar(p_stress, q_stress, exit, M, p0, beta, mu, K, xi, rma_prefac, particles.eps_pl_vol[p]);
+                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!");

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; p<Np; p++){
 
         //// Position is updated according to PIC velocities

src/simulation/save_data.cpp

@@ -25,12 +25,7 @@ void Simulation::saveParticleData(std::string extra){
 
         T pressure  = -tau.trace() / dim;
         TM tau_dev = tau + pressure * TM::Identity();
-
-        T devstress;
-        if (use_mises_q)
-            devstress = std::sqrt(3.0/2.0 * selfDoubleDot(tau_dev));
-        else
-            devstress = std::sqrt(1.0/2.0 * selfDoubleDot(tau_dev));
+        T devstress = q_prefac * std::sqrt(selfDoubleDot(tau_dev));
 
         pressure_vec[p] = pressure;
         devstress_vec[p] = devstress;

src/simulation/simulation.cpp

@@ -83,13 +83,16 @@ void Simulation::simulate(){
         return;
     }
 
-    #if DIMENSION == 3
+    if (dim == 3){
         debug("This is a 3D simulation.");
-    #elif DIMENSION == 2
+    }
+    else if (dim == 2){
         debug("This is a 2D simulation.");
-    #else
-        #error Unsupported spline degree
-    #endif
+    }
+    else{
+        debug("Unsupported spline degree");
+        return;
+    }
 
     #if SPLINEDEG == 3
       apicDinverse = 3.0/(dx*dx);
@@ -121,17 +124,9 @@ void Simulation::simulate(){
     one_over_dx = 1.0 / dx;
     one_over_dx_square = one_over_dx * one_over_dx;
 
-    if (use_mises_q){
-        q_prefac = std::sqrt(3.0)/std::sqrt(2.0);
-        d_prefac = std::sqrt(2.0)/std::sqrt(3.0);
-    } 
-    else {
-        q_prefac = 1.0 / std::sqrt(2.0);
-        d_prefac = std::sqrt(2.0);
-    }
-    e_mu_prefac = 2*q_prefac          * mu;
-    f_mu_prefac = 2*q_prefac/d_prefac * mu;
-    rma_prefac  = 2*q_prefac*q_prefac;
+    d_prefac = 1 / q_prefac;
+    e_mu_prefac = 2*q_prefac            * mu;
+    f_mu_prefac = 2*q_prefac * q_prefac * mu;
 
     fac_Q = I_ref / (grain_diameter*std::sqrt(rho_s));