From 5cee0e3e16300e311d9730235fee0aa6d27ff348 Mon Sep 17 00:00:00 2001 From: Lars Blatny Date: Thu, 14 Aug 2025 14:34:44 +0200 Subject: [PATCH 1/6] cleaned up prefactors, removed duplicate rma_prefac --- src/plasticity_helpers/mcc_rma_explicit.hpp | 12 +++---- .../mcc_rma_explicit_onevar.hpp | 3 +- .../mcc_rma_implicit_exponential.hpp | 12 +++---- .../mcc_rma_implicit_exponential_onevar.hpp | 3 +- .../mcc_rma_implicit_sinh_onevar.hpp | 3 +- src/simulation/plasticity.cpp | 32 +++++++++---------- src/simulation/simulation.cpp | 8 ++--- src/simulation/simulation.hpp | 1 - 8 files changed, 34 insertions(+), 40 deletions(-) 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 a4bb14a..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 */ { diff --git a/src/plasticity_helpers/mcc_rma_implicit_sinh_onevar.hpp b/src/plasticity_helpers/mcc_rma_implicit_sinh_onevar.hpp index 040f506..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 */ { diff --git a/src/simulation/plasticity.cpp b/src/simulation/plasticity.cpp index 1c46c16..b207156 100644 --- a/src/simulation/plasticity.cpp +++ b/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 = } 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!"); diff --git a/src/simulation/simulation.cpp b/src/simulation/simulation.cpp index c15f4a9..cf97765 100644 --- a/src/simulation/simulation.cpp +++ b/src/simulation/simulation.cpp @@ -123,15 +123,13 @@ void Simulation::simulate(){ 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)); diff --git a/src/simulation/simulation.hpp b/src/simulation/simulation.hpp index 25dc426..8923f18 100644 --- a/src/simulation/simulation.hpp +++ b/src/simulation/simulation.hpp @@ -177,7 +177,6 @@ class Simulation{ 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; From 4c640e7189c46a5fa3eb4dafdb61b0b80b1b1a6b Mon Sep 17 00:00:00 2001 From: Lars Blatny Date: Thu, 14 Aug 2025 14:44:27 +0200 Subject: [PATCH 2/6] stored extent of vdb sampled area in Lx, Ly and Lz --- README.md | 2 +- src/sampling/sampling_particles_vdb.hpp | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 0dc1425..33e74ee 100644 --- a/README.md +++ b/README.md @@ -210,7 +210,7 @@ This is a non-exhaustive list of parameters and options (of the `Simulation` cla | `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. diff --git a/src/sampling/sampling_particles_vdb.hpp b/src/sampling/sampling_particles_vdb.hpp index 2723b91..1973a6e 100644 --- a/src/sampling/sampling_particles_vdb.hpp +++ b/src/sampling/sampling_particles_vdb.hpp @@ -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,6 +52,9 @@ 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; + + sim.Lx = L(0); + sim.Ly = L(1); #endif // DIMENSION debug(" Number of square samples: ", square_samples.size()); From f27e2f817d6a453bc8d8ee551a4850f8bad86606 Mon Sep 17 00:00:00 2001 From: Lars Blatny Date: Thu, 21 Aug 2025 11:53:03 +0200 Subject: [PATCH 3/6] removed used of DIMENSION for non-omp reasons --- src/sampling/sampling_particles.hpp | 2 +- src/sampling/sampling_particles_vdb.hpp | 6 +++--- src/simulation/simulation.cpp | 13 ++++++++----- 3 files changed, 12 insertions(+), 9 deletions(-) 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 1973a6e..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..."); @@ -55,7 +55,7 @@ void sampleParticlesFromVdb(S& sim, ObjectVdb& obj, T kRadius, T ppc = 6) sim.Lx = L(0); sim.Ly = L(1); - #endif // DIMENSION + #endif debug(" Number of square samples: ", square_samples.size()); debug(" dx set to ", sim.dx); @@ -67,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/simulation.cpp b/src/simulation/simulation.cpp index cf97765..3da9c6e 100644 --- a/src/simulation/simulation.cpp +++ b/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); From 73977589daa69ee9da1ccdfe634404d2b745708d Mon Sep 17 00:00:00 2001 From: Lars Blatny Date: Thu, 21 Aug 2025 12:13:19 +0200 Subject: [PATCH 4/6] removed use_mises_q option, use instead directly the q_prefac variable --- README.md | 2 +- src/mpm.cpp | 2 +- src/simulation/save_data.cpp | 7 +------ src/simulation/simulation.cpp | 6 ------ src/simulation/simulation.hpp | 5 +++-- tests/tests.cpp | 15 ++++++--------- 6 files changed, 12 insertions(+), 25 deletions(-) diff --git a/README.md b/README.md index 33e74ee..7645b3f 100644 --- a/README.md +++ b/README.md @@ -208,7 +208,6 @@ 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(...)`. 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 @@ -216,6 +215,7 @@ This is a non-exhaustive list of parameters and options (of the `Simulation` cla | `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: diff --git a/src/mpm.cpp b/src/mpm.cpp index 84bf98e..8c753d2 100644 --- a/src/mpm.cpp +++ b/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 diff --git a/src/simulation/save_data.cpp b/src/simulation/save_data.cpp index 252dd92..dcb94c8 100644 --- a/src/simulation/save_data.cpp +++ b/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; diff --git a/src/simulation/simulation.cpp b/src/simulation/simulation.cpp index 3da9c6e..026da4a 100644 --- a/src/simulation/simulation.cpp +++ b/src/simulation/simulation.cpp @@ -124,12 +124,6 @@ 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); - } - else { - q_prefac = 1.0 / std::sqrt(2.0); - } d_prefac = 1 / q_prefac; e_mu_prefac = 2*q_prefac * mu; f_mu_prefac = 2*q_prefac * q_prefac * mu; diff --git a/src/simulation/simulation.hpp b/src/simulation/simulation.hpp index 8923f18..9a3acac 100644 --- a/src/simulation/simulation.hpp +++ b/src/simulation/simulation.hpp @@ -75,7 +75,6 @@ class Simulation{ PlasticModel plastic_model = PlasticModel::NoPlasticity; HardeningLaw hardening_law = HardeningLaw::ExpoImpl; bool use_pradhana = true; - bool use_mises_q = false; T E = 1e6; // Young's modulus (3D) T nu = 0.3; // Poisson's ratio (3D) T stress_tolerance = 1e-5; @@ -105,6 +104,9 @@ class Simulation{ T mu_1 = std::tan(20.9*M_PI/180.0); T mu_2 = std::tan(32.8*M_PI/180.0);; + // 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; @@ -173,7 +175,6 @@ class Simulation{ 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 diff --git a/tests/tests.cpp b/tests/tests.cpp index d05fc7c..07c3bcf 100644 --- a/tests/tests.cpp +++ b/tests/tests.cpp @@ -346,8 +346,7 @@ TEST(BoundaryTest, MIBF) { 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; @@ -398,8 +397,7 @@ TEST(BoundaryTest, MIBF) { 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; @@ -454,10 +452,11 @@ TEST(BoundaryTest, MIBF) { 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; @@ -641,7 +640,6 @@ TEST(CollapseTest, DruckerPragerOne) { sim.elastic_model = ElasticModel::Hencky; sim.plastic_model = PlasticModel::DPSoft; - sim.use_mises_q = false; sim.use_pradhana = true; sim.xi = 0; @@ -707,7 +705,6 @@ TEST(CollapseTest, DruckerPragerTwo) { sim.elastic_model = ElasticModel::Hencky; sim.plastic_model = PlasticModel::DPVisc; - sim.use_mises_q = false; sim.use_pradhana = true; sim.perzyna_visc = 0; From a35189bd58a4c060b9f4162c210d73fdeffd68e5 Mon Sep 17 00:00:00 2001 From: Lars Blatny Date: Thu, 21 Aug 2025 12:21:56 +0200 Subject: [PATCH 5/6] update example after latest change of mises_q --- examples/collapse.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/collapse.cpp b/examples/collapse.cpp index 84bf98e..8c753d2 100644 --- a/examples/collapse.cpp +++ b/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 From fcce80a33e2b30bca62097340e3f7219d317cea3 Mon Sep 17 00:00:00 2001 From: Lars Blatny Date: Mon, 25 Aug 2025 18:46:36 +0200 Subject: [PATCH 6/6] openmp nowaits added --- src/simulation/explicit_euler_update.cpp | 2 +- src/simulation/g2p.cpp | 2 +- src/simulation/musl.cpp | 4 ++-- src/simulation/p2g.cpp | 2 +- src/simulation/position_update.cpp | 2 +- 5 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/simulation/explicit_euler_update.cpp b/src/simulation/explicit_euler_update.cpp index 3e0fe27..d7b003a 100644 --- a/src/simulation/explicit_euler_update.cpp +++ b/src/simulation/explicit_euler_update.cpp @@ -16,7 +16,7 @@ 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]; 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/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; p