From 8c6da72141c90dd37cea8e9ac9795f889094d8e7 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Fri, 22 Dec 2023 09:25:14 +0100 Subject: [PATCH 01/24] Start on JHR model --- src/jhr.rs | 292 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 292 insertions(+) create mode 100644 src/jhr.rs diff --git a/src/jhr.rs b/src/jhr.rs new file mode 100644 index 0000000..cb740b3 --- /dev/null +++ b/src/jhr.rs @@ -0,0 +1,292 @@ + +use crate::interfaces::{ConstitutiveModel, QDim, QValueInput, QValueOutput, Q}; +use crate::stress_strain::{ + deviatoric, mandel_decomposition, mandel_rate_from_velocity_gradient, MANDEL_IDENTITY, +}; + +use nalgebra::{SMatrix, DVectorView}; +use std::cmp; +use std::collections::HashMap; + +#[derive(Debug)] +pub struct JHR3D { + _density: f64, + _a: f64, + _b: f64, + _n: f64, + _f_c: f64, + _s_max: f64, + _shear_modulus: f64, + _d_1: f64, + _d_2: f64, + _eps_f_min: f64, + _p_crush: f64, + _mu_crush: f64, + _p_lock: f64, + _mu_lock: f64, + _k_1: f64, + _k_2: f64, + _k_3: f64, + _t: f64, +} + +impl ConstitutiveModel for JHR3D { + fn new(parameters: &HashMap) -> Option { + Some(Self { + _density: *parameters.get("density")?, + _a: *parameters.get("a")?, + _b: *parameters.get("b")?, + _n: *parameters.get("n")?, + _f_c: *parameters.get("f_c")?, + _s_max: *parameters.get("s_max")?, + _shear_modulus: *parameters.get("shear_modulus")?, + _d_1: *parameters.get("d_1")?, + _d_2: *parameters.get("d_2")?, + _eps_f_min: *parameters.get("eps_f_min")?, + _p_crush: *parameters.get("p_crush")?, + _mu_crush: *parameters.get("mu_crush")?, + _p_lock: *parameters.get("p_lock")?, + _mu_lock: *parameters.get("mu_lock")?, + _k_1: *parameters.get("k_1")?, + _k_2: *parameters.get("k_2")?, + _k_3: *parameters.get("k_3")?, + _t: *parameters.get("t")?, + }) + } + fn evaluate_ip(&self, ip: usize, del_t: f64, input: &QValueInput, output: &mut QValueOutput) { + let velocity_gradient = input + .get_tensor::<{ Q::VelocityGradient.dim() }, { Q::VelocityGradient.size() }>( + Q::VelocityGradient, + ip, + ); + let d_eps = mandel_rate_from_velocity_gradient(&velocity_gradient); + let (mut d_eps_vol, d_eps_dev) = mandel_decomposition(&d_eps); + d_eps_vol *= -1.0; + + let sigma_0 = input.get_vector::<{ Q::MandelStress.size() }>(Q::MandelStress, ip); + + let mut del_lambda = 0.0; + + let damage_0 = input.get_scalar(Q::Damage, ip); + let mut damage_1 = damage_0; + + let (p_0, s_0) = mandel_decomposition(&sigma_0); + let p_0 = p_0 - input.get_scalar(Q::BulkViscosity, ip); + + let s_tr = s_0 + 2. * self._shear_modulus * d_eps_dev * del_t; + let s_tr_eq = (1.5 * s_tr.norm_squared()).sqrt(); + let d_eps_eq = ((2. / 3.) * d_eps.norm_squared()).sqrt(); + let mut alpha; + + let p_s = p_0 / self.parameters.PHEL; + let t_s = self.parameters.T / self.parameters.PHEL; + let mut rate_factor = 1.; + + let fracture_surface = + (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) + .max(0.0); + let residual_surface = + (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); + if d_eps_eq >= self.parameters.EPS0 { + rate_factor += self.parameters.C * (d_eps_eq / self.parameters.EPS0).ln(); + } + let yield_surface = rate_factor * { + if damage_0 == 0.0 { + fracture_surface + } else { + fracture_surface * (1. - damage_0) + damage_0 * residual_surface + } + }; + if s_tr_eq > yield_surface { + let e_p_f = (self.parameters.D1 * (p_s + t_s).powf(self.parameters.D2)) + .max(self.parameters.EFMIN); + + del_lambda = (s_tr_eq - yield_surface) / (3. * self.parameters.SHEAR_MODULUS); + alpha = yield_surface / s_tr_eq; + + damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); + output.set_scalar(Q::Damage, ip, damage_1); + } else { + alpha = 1.0; + } + + // /*********************************************************************** + // * UPDATE DENSITY + // * The density is updated using the explicit midpoint rule for the + // * deformation gradient. + // TODO: Move this since, it will be calculated outside of the constitutive model + // **********************************************************************/ + let f1 = del_t / 2. * 3. * d_eps_vol; + let density_0 = input.get_scalar(Q::Density, ip); + let density_1 = density_0 * (1. - f1) / (1. + f1); + debug_assert!( + density_1 > 0.0, + "Negative density encountered in JH2 model: {}", + density_1 + ); + output.set_scalar(Q::Density, ip, density_1); + + let mu = density_1 / self.parameters.RHO - 1.; + + let p_1 = { + if mu > 0.0 { + self.parameters.K1 * mu + + self.parameters.K2 * mu.powi(2) + + self.parameters.K3 * mu.powi(3) + + input.get_scalar(Q::BulkingPressure, ip) + } else { + let p_trial = self.parameters.K1 * mu; + let p_damaged = -self.parameters.T * (1. - damage_1); + if p_trial > p_damaged { + p_trial + } else { + p_damaged + } + } + }; + if damage_1 > damage_0 { + let y_old = damage_0 * residual_surface + (1. - damage_0) * fracture_surface; + let y_new = damage_1 * residual_surface + (1. - damage_1) * fracture_surface; + let u_old = y_old.powi(2) / (6. * self.parameters.SHEAR_MODULUS); + let u_new = y_new.powi(2) / (6. * self.parameters.SHEAR_MODULUS); + + let del_u = u_old - u_new; + + let del_p_0 = input.get_scalar(Q::BulkingPressure, ip); + let del_p = -self.parameters.K1 * mu + + ((self.parameters.K1 * mu + del_p_0).powi(2) + + 2. * self.parameters.BETA * self.parameters.K1 * del_u) + .sqrt(); + output.set_scalar(Q::BulkingPressure, ip, del_p); + } + + // Calculate bulk viscosity + let l = input.get_scalar(Q::CellDiameter, ip); + let tr_deps = d_eps_vol * 3.; + let c = (self.parameters.K1 / self.parameters.RHO).sqrt(); + let q_1 = { + if tr_deps < 0.0 { + l * density_1 * (1.5 * l * tr_deps.powi(2) + c * 0.06 * tr_deps.abs()) + } else { + 0.0 + } + }; + output.set_scalar(Q::BulkViscosity, ip, q_1); + // /*********************************************************************** + // * Combine deviatoric and volumetric stresses + // **********************************************************************/ + let s_1 = s_tr * alpha; + let sigma_1 = s_1 - MANDEL_IDENTITY * (p_1 + q_1); + output.set_vector(Q::MandelStress, ip, sigma_1); + + // *********************************************************************** + // Update optional output variables if needed + // ********************************************************************** + + if output.is_some(Q::EqStrainRate) { + output.set_scalar(Q::EqStrainRate, ip, d_eps_eq); + } + if output.is_some(Q::MandelStrainRate) { + output.set_vector(Q::MandelStrainRate, ip, d_eps); + } + if output.is_some(Q::MisesStress) { + output.set_scalar(Q::MisesStress, ip, alpha * s_tr_eq); + } + if output.is_some(Q::Pressure) { + output.set_scalar(Q::Pressure, ip, p_1); + } + + let elastic_rate = + - (1. - alpha) / (2. * self.parameters.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; + + let density_mid = 0.5 * (density_0 + density_1); + if output.is_some(Q::InternalPlasticEnergy) && input.is_some(Q::InternalPlasticEnergy) { + let s_mid = 0.5 * (s_0 + s_1); + let deviatoric_rate = d_eps_dev - elastic_rate; + let e_0 = input.get_scalar(Q::InternalPlasticEnergy, ip); + let e_1 = e_0 + + del_t / density_mid * (s_mid.dot(&deviatoric_rate)); + output.set_scalar(Q::InternalPlasticEnergy, ip, e_1); + } + if output.is_some(Q::InternalElasticEnergy) && input.is_some(Q::InternalElasticEnergy) { + let s_mid = 0.5 * (s_0 + s_1); + let p_mid = -0.5 * (p_0 + p_1); + let deviatoric_rate = elastic_rate; + let e_0 = input.get_scalar(Q::InternalElasticEnergy, ip); + let e_1 = e_0 + + del_t / density_mid + * (s_mid.dot(&deviatoric_rate) + 3. * d_eps_vol * p_mid); + output.set_scalar(Q::InternalElasticEnergy, ip, e_1); + } + if output.is_some(Q::InternalEnergy) && input.is_some(Q::InternalEnergy) { + let e_0 = input.get_scalar(Q::InternalEnergy, ip); + let sigma_mid = 0.5 * (sigma_0 + sigma_1); + let e_1 = e_0 + del_t / density_mid * sigma_mid.dot(&d_eps); + output.set_scalar(Q::InternalEnergy, ip, e_1); + } + if output.is_some(Q::InternalHeatingEnergy) && input.is_some(Q::InternalHeatingEnergy) { + let e_0 = input.get_scalar(Q::InternalHeatingEnergy, ip); + let q_mid = - 0.5 * (q_1 + input.get_scalar(Q::BulkViscosity, ip)); + let e_1 = e_0 + del_t / density_mid * 3. * q_mid * d_eps_vol; + output.set_scalar(Q::InternalHeatingEnergy, ip, e_1); + } + if output.is_some(Q::EqPlasticStrain) && input.is_some(Q::EqPlasticStrain) { + output.set_scalar( + Q::EqPlasticStrain, + ip, + input.get_scalar(Q::EqPlasticStrain, ip) + del_lambda, + ); + } + } + + /// Returns the physical quantities that are required as input for the + /// constitutive model together with their dimensions. + fn define_input(&self) -> HashMap { + HashMap::from([ + (Q::VelocityGradient, QDim::SquareTensor(3)), + (Q::CellDiameter, QDim::Scalar), + ]) + } + + /// Returns the physical quantities that are needed as internal variables + /// for the constitutive model together with their dimensions. These Variables are + /// stored both in in the input and the output. + fn define_history(&self) -> HashMap { + HashMap::from([ + (Q::MandelStress, QDim::Vector(6)), + (Q::Damage, QDim::Scalar), + (Q::BulkingPressure, QDim::Scalar), + (Q::Density, QDim::Scalar), + (Q::BulkViscosity, QDim::Scalar), + ]) + } + + /// Returns the physical quantities that are needed as output, but are not + /// necessarily needed in oredr to calculate the constitutive model. An example is + /// the consistent tangent which is not needed for the calculation of the stresses + /// and is therefore purely an output quantity. + fn define_output(&self) -> HashMap { + HashMap::from([(Q::MandelStress, QDim::Vector(6))]) + } + + /// Returns the physical quantities that are optional output of the constitutive + /// model. These quantities are not needed for the calculation of the stresses + /// but can be useful for postprocessing. + fn define_optional_output(&self) -> HashMap { + HashMap::from([ + (Q::EqStrainRate, QDim::Scalar), + (Q::MandelStrainRate, QDim::Vector(6)), + (Q::MisesStress, QDim::Scalar), + (Q::Pressure, QDim::Scalar), + ]) + } + fn define_optional_history(&self) -> HashMap { + HashMap::from([ + (Q::EqPlasticStrain, QDim::Scalar), + (Q::InternalPlasticEnergy, QDim::Scalar), + (Q::InternalElasticEnergy, QDim::Scalar), + (Q::InternalEnergy, QDim::Scalar), + (Q::InternalHeatingEnergy, QDim::Scalar), + ]) + } +} From 01204459322018d09b99fadf6dbb578f525e371e Mon Sep 17 00:00:00 2001 From: srosenbu Date: Fri, 22 Dec 2023 12:23:10 +0100 Subject: [PATCH 02/24] Parameter for tensile collapse of yield surface --- src/gradient_jh2.rs | 6 +- src/jh2.rs | 6 +- src/jhr.rs | 130 ++++++++++++++++++++------------------------ src/lib.rs | 2 + 4 files changed, 68 insertions(+), 76 deletions(-) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index d61e8f6..872f9e0 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -35,6 +35,7 @@ impl ConstitutiveModel for GradientJH23D { BETA: *parameters.get("BETA").unwrap(), EFMIN: *parameters.get("EFMIN").unwrap(), DMAX: *parameters.get("DMAX").unwrap_or(&1.0), + REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, }) } @@ -72,9 +73,10 @@ impl ConstitutiveModel for GradientJH23D { let damage_0 = input.get_scalar(Q::Damage, ip); let damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); output.set_scalar(Q::Damage, ip, damage_1); - + + let t_s_factor = (1.-damage_1).powf(self.parameters.REDUCE_T); let fracture_surface = - (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) + (self.parameters.A * (p_s + t_s * t_s_factor).powf(self.parameters.N) * self.parameters.SIGMAHEL) .max(0.0); let residual_surface = (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); diff --git a/src/jh2.rs b/src/jh2.rs index f0498e0..f95e5fa 100644 --- a/src/jh2.rs +++ b/src/jh2.rs @@ -25,6 +25,7 @@ pub struct JH2ConstParameters { pub BETA: f64, pub EFMIN: f64, pub DMAX: f64, + pub REDUCE_T: f64, } #[derive(Debug)] pub struct JH23D { @@ -54,6 +55,7 @@ impl ConstitutiveModel for JH23D { BETA: *parameters.get("BETA").unwrap(), EFMIN: *parameters.get("EFMIN").unwrap(), DMAX: *parameters.get("DMAX").unwrap_or(&1.0), + REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, }) } @@ -85,9 +87,9 @@ impl ConstitutiveModel for JH23D { let p_s = p_0 / self.parameters.PHEL; let t_s = self.parameters.T / self.parameters.PHEL; let mut rate_factor = 1.; - + let t_s_factor = (1.-damage_0).powf(self.parameters.REDUCE_T); let fracture_surface = - (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) + (self.parameters.A * (p_s + t_s * t_s_factor).powf(self.parameters.N) * self.parameters.SIGMAHEL) .max(0.0); let residual_surface = (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); diff --git a/src/jhr.rs b/src/jhr.rs index cb740b3..b96c1f6 100644 --- a/src/jhr.rs +++ b/src/jhr.rs @@ -10,47 +10,49 @@ use std::collections::HashMap; #[derive(Debug)] pub struct JHR3D { - _density: f64, - _a: f64, - _b: f64, - _n: f64, - _f_c: f64, - _s_max: f64, - _shear_modulus: f64, - _d_1: f64, - _d_2: f64, - _eps_f_min: f64, - _p_crush: f64, - _mu_crush: f64, - _p_lock: f64, - _mu_lock: f64, - _k_1: f64, - _k_2: f64, - _k_3: f64, - _t: f64, + RHO: f64, + SHEAR_MODULUS: f64, + A: f64, + B: f64, + C: f64, + M: f64, + N: f64, + EPS0: f64, + T: f64, + SIGMAHEL: f64, + PHEL: f64, + D1: f64, + D2: f64, + K1: f64, + K2: f64, + K3: f64, + BETA: f64, + EFMIN: f64, + DMAX: f64, } impl ConstitutiveModel for JHR3D { fn new(parameters: &HashMap) -> Option { Some(Self { - _density: *parameters.get("density")?, - _a: *parameters.get("a")?, - _b: *parameters.get("b")?, - _n: *parameters.get("n")?, - _f_c: *parameters.get("f_c")?, - _s_max: *parameters.get("s_max")?, - _shear_modulus: *parameters.get("shear_modulus")?, - _d_1: *parameters.get("d_1")?, - _d_2: *parameters.get("d_2")?, - _eps_f_min: *parameters.get("eps_f_min")?, - _p_crush: *parameters.get("p_crush")?, - _mu_crush: *parameters.get("mu_crush")?, - _p_lock: *parameters.get("p_lock")?, - _mu_lock: *parameters.get("mu_lock")?, - _k_1: *parameters.get("k_1")?, - _k_2: *parameters.get("k_2")?, - _k_3: *parameters.get("k_3")?, - _t: *parameters.get("t")?, + RHO: *parameters.get("RHO")?, + SHEAR_MODULUS: *parameters.get("SHEAR_MODULUS")?, + A: *parameters.get("A")?, + B: *parameters.get("B")?, + C: *parameters.get("C")?, + M: *parameters.get("M")?, + N: *parameters.get("N").unwrap(), + EPS0: *parameters.get("EPS0").unwrap(), + T: *parameters.get("T").unwrap(), + SIGMAHEL: *parameters.get("SIGMAHEL").unwrap(), + PHEL: *parameters.get("PHEL").unwrap(), + D1: *parameters.get("D1").unwrap(), + D2: *parameters.get("D2").unwrap(), + K1: *parameters.get("K1").unwrap(), + K2: *parameters.get("K2").unwrap(), + K3: *parameters.get("K3").unwrap(), + BETA: *parameters.get("BETA").unwrap(), + EFMIN: *parameters.get("EFMIN").unwrap(), + DMAX: *parameters.get("DMAX").unwrap_or(&1.0), }) } fn evaluate_ip(&self, ip: usize, del_t: f64, input: &QValueInput, output: &mut QValueOutput) { @@ -73,22 +75,22 @@ impl ConstitutiveModel for JHR3D { let (p_0, s_0) = mandel_decomposition(&sigma_0); let p_0 = p_0 - input.get_scalar(Q::BulkViscosity, ip); - let s_tr = s_0 + 2. * self._shear_modulus * d_eps_dev * del_t; + let s_tr = s_0 + 2. * self.SHEAR_MODULUS * d_eps_dev * del_t; let s_tr_eq = (1.5 * s_tr.norm_squared()).sqrt(); let d_eps_eq = ((2. / 3.) * d_eps.norm_squared()).sqrt(); let mut alpha; - let p_s = p_0 / self.parameters.PHEL; - let t_s = self.parameters.T / self.parameters.PHEL; + let p_s = p_0 / self.PHEL; + let t_s = self.T / self.PHEL; let mut rate_factor = 1.; let fracture_surface = - (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) + (self.A * (p_s + t_s*(1.-damage_0)).powf(self.N) * self.SIGMAHEL) .max(0.0); let residual_surface = - (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); - if d_eps_eq >= self.parameters.EPS0 { - rate_factor += self.parameters.C * (d_eps_eq / self.parameters.EPS0).ln(); + (self.B * (p_s).powf(self.M) * self.SIGMAHEL).max(0.0); + if d_eps_eq >= self.EPS0 { + rate_factor += self.C * (d_eps_eq / self.EPS0).ln(); } let yield_surface = rate_factor * { if damage_0 == 0.0 { @@ -98,13 +100,13 @@ impl ConstitutiveModel for JHR3D { } }; if s_tr_eq > yield_surface { - let e_p_f = (self.parameters.D1 * (p_s + t_s).powf(self.parameters.D2)) - .max(self.parameters.EFMIN); + let e_p_f = (self.D1 * (p_s + t_s).powf(self.D2)) + .max(self.EFMIN); - del_lambda = (s_tr_eq - yield_surface) / (3. * self.parameters.SHEAR_MODULUS); + del_lambda = (s_tr_eq - yield_surface) / (3. * self.SHEAR_MODULUS); alpha = yield_surface / s_tr_eq; - damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); + damage_1 = (damage_0 + del_lambda / e_p_f).min(self.DMAX); output.set_scalar(Q::Damage, ip, damage_1); } else { alpha = 1.0; @@ -126,17 +128,17 @@ impl ConstitutiveModel for JHR3D { ); output.set_scalar(Q::Density, ip, density_1); - let mu = density_1 / self.parameters.RHO - 1.; + let mu = density_1 / self.RHO - 1.; let p_1 = { if mu > 0.0 { - self.parameters.K1 * mu - + self.parameters.K2 * mu.powi(2) - + self.parameters.K3 * mu.powi(3) - + input.get_scalar(Q::BulkingPressure, ip) + self.K1 * mu + + self.K2 * mu.powi(2) + + self.K3 * mu.powi(3) + //+ input.get_scalar(Q::BulkingPressure, ip) } else { - let p_trial = self.parameters.K1 * mu; - let p_damaged = -self.parameters.T * (1. - damage_1); + let p_trial = self.K1 * mu; + let p_damaged = -self.T * (1. - damage_1); if p_trial > p_damaged { p_trial } else { @@ -144,26 +146,11 @@ impl ConstitutiveModel for JHR3D { } } }; - if damage_1 > damage_0 { - let y_old = damage_0 * residual_surface + (1. - damage_0) * fracture_surface; - let y_new = damage_1 * residual_surface + (1. - damage_1) * fracture_surface; - let u_old = y_old.powi(2) / (6. * self.parameters.SHEAR_MODULUS); - let u_new = y_new.powi(2) / (6. * self.parameters.SHEAR_MODULUS); - - let del_u = u_old - u_new; - - let del_p_0 = input.get_scalar(Q::BulkingPressure, ip); - let del_p = -self.parameters.K1 * mu - + ((self.parameters.K1 * mu + del_p_0).powi(2) - + 2. * self.parameters.BETA * self.parameters.K1 * del_u) - .sqrt(); - output.set_scalar(Q::BulkingPressure, ip, del_p); - } // Calculate bulk viscosity let l = input.get_scalar(Q::CellDiameter, ip); let tr_deps = d_eps_vol * 3.; - let c = (self.parameters.K1 / self.parameters.RHO).sqrt(); + let c = (self.K1 / self.RHO).sqrt(); let q_1 = { if tr_deps < 0.0 { l * density_1 * (1.5 * l * tr_deps.powi(2) + c * 0.06 * tr_deps.abs()) @@ -197,7 +184,7 @@ impl ConstitutiveModel for JHR3D { } let elastic_rate = - - (1. - alpha) / (2. * self.parameters.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; + - (1. - alpha) / (2. * self.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; let density_mid = 0.5 * (density_0 + density_1); if output.is_some(Q::InternalPlasticEnergy) && input.is_some(Q::InternalPlasticEnergy) { @@ -255,7 +242,6 @@ impl ConstitutiveModel for JHR3D { HashMap::from([ (Q::MandelStress, QDim::Vector(6)), (Q::Damage, QDim::Scalar), - (Q::BulkingPressure, QDim::Scalar), (Q::Density, QDim::Scalar), (Q::BulkViscosity, QDim::Scalar), ]) diff --git a/src/lib.rs b/src/lib.rs index 3a68320..8378f35 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -3,6 +3,7 @@ use std::collections::HashMap; use crate::interfaces::{ConstitutiveModel, QDim, QValueInput, QValueOutput, Q}; use crate::jh2::JH23D; +use crate::jhr::JHR3D; //use crate::jh_concrete::JHConcrete3D; use crate::generic_jh2::GenericJH23D; use crate::gradient_jh2::GradientJH23D; @@ -19,6 +20,7 @@ pub mod jh2; //pub mod jh_concrete; pub mod generic_jh2; pub mod gradient_jh2; +pub mod jhr; pub mod smallstrain; pub mod stress_strain; From 3840e7276ce0171d25ced863af83440f147a5f25 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Fri, 26 Jan 2024 16:04:41 +0100 Subject: [PATCH 03/24] Functional programming --- src/jh2.rs | 1 + src/lib.rs | 19 +++++- src/smallstrain.rs | 140 +++++++++++++++++++++++++++++++++------------ 3 files changed, 121 insertions(+), 39 deletions(-) diff --git a/src/jh2.rs b/src/jh2.rs index f95e5fa..219f514 100644 --- a/src/jh2.rs +++ b/src/jh2.rs @@ -87,6 +87,7 @@ impl ConstitutiveModel for JH23D { let p_s = p_0 / self.parameters.PHEL; let t_s = self.parameters.T / self.parameters.PHEL; let mut rate_factor = 1.; + let t_s_factor = (1.-damage_0).powf(self.parameters.REDUCE_T); let fracture_surface = (self.parameters.A * (p_s + t_s * t_s_factor).powf(self.parameters.N) * self.parameters.SIGMAHEL) diff --git a/src/lib.rs b/src/lib.rs index 8378f35..2c6d6a1 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -8,6 +8,7 @@ use crate::jhr::JHR3D; use crate::generic_jh2::GenericJH23D; use crate::gradient_jh2::GradientJH23D; use crate::smallstrain::linear_elastic::LinearElastic3D; +use crate::smallstrain::{evaluate_model, elasticity_3d}; //use crate::stress_strain; use nalgebra::{Const, DVectorView, DVectorViewMut, Dyn, SMatrix}; use numpy::{PyReadonlyArray1, PyReadwriteArray1}; @@ -393,7 +394,23 @@ fn py_jaumann_rotation_expensive( stress_strain::jaumann_rotation_expensive(del_t, &velocity_gradient, &mut stress); Ok(()) } - +#[pyfunction(name="evaluate_elasticity_3d")] +fn py_evaluate_elasticity_3d( + del_t: f64, + stress: PyReadwriteArray1, + del_strain: PyReadonlyArray1, + parameters: PyReadonlyArray1, + history: PyReadwriteArray1, + tangent: PyReadwriteArray1, +) -> PyResult<()> { + let stress = stress.as_slice_mut()?; + let del_strain = del_strain.as_slice()?; + let parameters = parameters.as_slice()?; + let history = history.as_slice_mut()?; + let mut tangent = tangent.as_slice_mut()?; + evaluate_model::<6,0,2>(&elasticity_3d, del_t, stress, del_strain, parameters, history, tangent); + Ok(()) +} #[pymodule] fn comfe(_py: Python, m: &PyModule) -> PyResult<()> { m.add_class::()?; diff --git a/src/smallstrain.rs b/src/smallstrain.rs index 8d65904..1561ed5 100644 --- a/src/smallstrain.rs +++ b/src/smallstrain.rs @@ -1,44 +1,108 @@ -use nalgebra::{DVectorView, DVectorViewMut}; +use nalgebra::{DVectorView, DVectorViewMut, SMatrix, SVector}; +use std::convert::TryInto; pub mod linear_elastic; -pub trait SmallStrainModel { - //fn new(parameters: &HashMap) -> Self; - fn evaluate_ip( - &self, - ip: usize, - del_t: f64, - stress: &mut DVectorViewMut, - del_strain: &DVectorView, - tangent: &mut DVectorViewMut, +// define a type for a specific function signature +pub type SmallStrainModelFn< + const STRESS_STRAIN_DIM: usize, + const HISTORY_DIM: usize, + const PARAMETER_DIM: usize, +> = fn( + del_t: f64, + stress: [f64; STRESS_STRAIN_DIM], + del_strain: [f64; STRESS_STRAIN_DIM], + parameters: [f64; PARAMETER_DIM], + history: [f64; HISTORY_DIM], +) -> ( + [f64; STRESS_STRAIN_DIM], + [[f64; STRESS_STRAIN_DIM]; STRESS_STRAIN_DIM], + [f64; HISTORY_DIM], +); + +pub fn elasticity_3d( + del_t: f64, + stress: [f64; 6], + del_strain: [f64; 6], + parameters: [f64; 2], + history: [f64; 0], +) -> ([f64; 6], [[f64; 6]; 6], [f64; 0]) { + let E = parameters[0]; + let nu = parameters[1]; + let mu = E / (2.0 * (1.0 + nu)); + let lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); + let c1 = lambda + 2.0 * mu; + let c2 = 2. * mu; + let D = SMatrix::::from([ + [c1, lambda, lambda, 0.0, 0.0, 0.0], + [lambda, c1, lambda, 0.0, 0.0, 0.0], + [lambda, lambda, c1, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, c2, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, c2, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, c2], + ]); + let stress_vector = SVector::::from(stress); + let del_strain_vector = SVector::::from(del_strain); + let new_stress = stress_vector + D * del_strain_vector; + let D_out = D.data.0; + let new_stress = new_stress.data.0[0]; + let history = history; + (new_stress, D_out, history) +} +pub struct SmallStrainModel< + const STRESS_STRAIN_DIM: usize, + const HISTORY_DIM: usize, + const PARAMETER_DIM: usize, +> { + pub model: SmallStrainModelFn, + pub parameters: [f64; PARAMETER_DIM], +} + +pub fn evaluate_model< + const STRESS_STRAIN_DIM: usize, + const HISTORY_DIM: usize, + const PARAMETER_DIM: usize, +>( + model: &SmallStrainModelFn, + del_t: f64, + stress: &mut [f64], + del_strain: &[f64], + parameters: &[f64], + history: &mut [f64], + tangent: &mut [f64], +) { + let parameters: [f64; PARAMETER_DIM] = parameters + .try_into() + .expect("Slice length does not match array length"); + + let stress_chunks = stress.chunks_exact_mut(STRESS_STRAIN_DIM); + let del_strain_chunks = del_strain.chunks_exact(STRESS_STRAIN_DIM); + let history_chunks = history.chunks_exact_mut(HISTORY_DIM); + let tangent_chunks = tangent.chunks_exact_mut(STRESS_STRAIN_DIM.pow(2)); + + assert!( + stress_chunks.len() == del_strain_chunks.len() + && stress_chunks.len() == history_chunks.len() + && stress_chunks.len() == tangent_chunks.len() ); - fn evaluate( - &self, - del_t: f64, - stress: &mut DVectorViewMut, - del_strain: &DVectorView, - tangent: &mut DVectorViewMut, - ) { - assert_eq!(stress.nrows(), del_strain.nrows()); - let n: usize = stress.nrows() / 6; - assert_eq!(stress.nrows(), n * 6); - for ip in 0..n { - self.evaluate_ip(ip, del_t, stress, del_strain, tangent); - } - } - fn evaluate_some( - &self, - del_t: f64, - stress: &mut DVectorViewMut, - del_strain: &DVectorView, - tangent: &mut DVectorViewMut, - ips: &[usize], - ) { - assert_eq!(stress.nrows(), del_strain.nrows()); - let n: usize = stress.nrows() / 6; - assert_eq!(stress.nrows(), n * 6); - for ip in ips { - self.evaluate_ip(*ip, del_t, stress, del_strain, tangent); - } + for (((stress_chunk, del_strain_chunk), history_chunk), tangent_chunk) in stress_chunks + .zip(del_strain_chunks) + .zip(history_chunks) + .zip(tangent_chunks) + { + let stress: [f64; STRESS_STRAIN_DIM] = stress_chunk + .try_into() + .expect("Slice length does not match array length"); + let del_strain: [f64; STRESS_STRAIN_DIM] = del_strain_chunk + .try_into() + .expect("Slice length does not match array length"); + let history: [f64; HISTORY_DIM] = history_chunk + .try_into() + .expect("Slice length does not match array length"); + let (new_stress, D, new_history) = model(del_t, stress, del_strain, parameters, history); + stress_chunk.copy_from_slice(&new_stress); + history_chunk.copy_from_slice(&new_history); + tangent_chunk.copy_from_slice(unsafe{std::slice::from_raw_parts(D.as_ptr() as *const f64, STRESS_STRAIN_DIM.pow(2))}); } + } From 1fdf0698b4b1d2a2d16c755881b27434ab4cb192 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 14 Feb 2024 14:10:31 +0100 Subject: [PATCH 04/24] exponential damage law --- src/gradient_jh2.rs | 17 +++++++++++++---- src/jh2.rs | 19 ++++++++++++++----- src/lib.rs | 34 +++++++++++++++++----------------- 3 files changed, 44 insertions(+), 26 deletions(-) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index 872f9e0..a02efdb 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -35,7 +35,8 @@ impl ConstitutiveModel for GradientJH23D { BETA: *parameters.get("BETA").unwrap(), EFMIN: *parameters.get("EFMIN").unwrap(), DMAX: *parameters.get("DMAX").unwrap_or(&1.0), - REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), + E_F: *parameters.get("E_F").unwrap_or(&0.0), + //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, }) } @@ -71,12 +72,20 @@ impl ConstitutiveModel for GradientJH23D { .max(self.parameters.EFMIN); let damage_0 = input.get_scalar(Q::Damage, ip); - let damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); + let mut damage_1 = 0.0; + if self.parameters.E_F > 0.0 { + let lambda_old = - self.parameters.E_F * (1. - damage_0).ln(); + let lambda_new = lambda_old + del_lambda_nonlocal; + damage_1 = 1. - (-lambda_new / self.parameters.E_F).exp(); + } else { + damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); + } + //let damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); output.set_scalar(Q::Damage, ip, damage_1); - let t_s_factor = (1.-damage_1).powf(self.parameters.REDUCE_T); + //let t_s_factor = (1.-damage_1).powf(self.parameters.REDUCE_T); let fracture_surface = - (self.parameters.A * (p_s + t_s * t_s_factor).powf(self.parameters.N) * self.parameters.SIGMAHEL) + (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) .max(0.0); let residual_surface = (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); diff --git a/src/jh2.rs b/src/jh2.rs index 219f514..22dfd8e 100644 --- a/src/jh2.rs +++ b/src/jh2.rs @@ -25,7 +25,8 @@ pub struct JH2ConstParameters { pub BETA: f64, pub EFMIN: f64, pub DMAX: f64, - pub REDUCE_T: f64, + pub E_F: f64, + //pub REDUCE_T: f64, } #[derive(Debug)] pub struct JH23D { @@ -55,7 +56,8 @@ impl ConstitutiveModel for JH23D { BETA: *parameters.get("BETA").unwrap(), EFMIN: *parameters.get("EFMIN").unwrap(), DMAX: *parameters.get("DMAX").unwrap_or(&1.0), - REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), + E_F: *parameters.get("E_F").unwrap_or(&0.0), + //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, }) } @@ -88,9 +90,9 @@ impl ConstitutiveModel for JH23D { let t_s = self.parameters.T / self.parameters.PHEL; let mut rate_factor = 1.; - let t_s_factor = (1.-damage_0).powf(self.parameters.REDUCE_T); + //let t_s_factor = (1.-damage_0).powf(self.parameters.REDUCE_T); let fracture_surface = - (self.parameters.A * (p_s + t_s * t_s_factor).powf(self.parameters.N) * self.parameters.SIGMAHEL) + (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) .max(0.0); let residual_surface = (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); @@ -111,7 +113,14 @@ impl ConstitutiveModel for JH23D { del_lambda = (s_tr_eq - yield_surface) / (3. * self.parameters.SHEAR_MODULUS); alpha = yield_surface / s_tr_eq; - damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); + if self.parameters.E_F > 0.0 { + let lambda_old = - self.parameters.E_F * (1. - damage_0).ln(); + let lambda_new = lambda_old + del_lambda; + damage_1 = 1. - (-lambda_new / self.parameters.E_F).exp(); + } else { + damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); + } + //damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); output.set_scalar(Q::Damage, ip, damage_1); } else { alpha = 1.0; diff --git a/src/lib.rs b/src/lib.rs index 2c6d6a1..eb1e74a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -394,23 +394,23 @@ fn py_jaumann_rotation_expensive( stress_strain::jaumann_rotation_expensive(del_t, &velocity_gradient, &mut stress); Ok(()) } -#[pyfunction(name="evaluate_elasticity_3d")] -fn py_evaluate_elasticity_3d( - del_t: f64, - stress: PyReadwriteArray1, - del_strain: PyReadonlyArray1, - parameters: PyReadonlyArray1, - history: PyReadwriteArray1, - tangent: PyReadwriteArray1, -) -> PyResult<()> { - let stress = stress.as_slice_mut()?; - let del_strain = del_strain.as_slice()?; - let parameters = parameters.as_slice()?; - let history = history.as_slice_mut()?; - let mut tangent = tangent.as_slice_mut()?; - evaluate_model::<6,0,2>(&elasticity_3d, del_t, stress, del_strain, parameters, history, tangent); - Ok(()) -} +// #[pyfunction(name="evaluate_elasticity_3d")] +// fn py_evaluate_elasticity_3d( +// del_t: f64, +// stress: PyReadwriteArray1, +// del_strain: PyReadonlyArray1, +// parameters: PyReadonlyArray1, +// history: PyReadwriteArray1, +// tangent: PyReadwriteArray1, +// ) -> PyResult<()> { +// let stress = stress.as_slice_mut()?; +// let del_strain = del_strain.as_slice()?; +// let parameters = parameters.as_slice()?; +// let history = history.as_slice_mut()?; +// let mut tangent = tangent.as_slice_mut()?; +// evaluate_model::<6,0,2>(&elasticity_3d, del_t, stress, del_strain, parameters, history, tangent); +// Ok(()) +// } #[pymodule] fn comfe(_py: Python, m: &PyModule) -> PyResult<()> { m.add_class::()?; From a15bf25075d8d8a2574232f175141c51a1945b84 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 14 Feb 2024 14:52:08 +0100 Subject: [PATCH 05/24] gradient model in terms of total nonlocal strain --- python/comfe/cdm.py | 15 +++++++-------- src/gradient_jh2.rs | 15 ++++++++++----- src/interfaces.rs | 3 +++ 3 files changed, 20 insertions(+), 13 deletions(-) diff --git a/python/comfe/cdm.py b/python/comfe/cdm.py index 749db35..5a1c952 100644 --- a/python/comfe/cdm.py +++ b/python/comfe/cdm.py @@ -361,7 +361,7 @@ class CDMNonlocal(NonlocalInterface): M: df.fem.Function # fields: dict[str, df.fem.Function] # q_fields: dict[str, df.fem.Function | ufl.core.expr.Expr] - rate_evaluator: QuadratureEvaluator + strain_evaluator: QuadratureEvaluator parameters: dict[str, float] form: df.fem.FormMetaClass t: float @@ -379,7 +379,7 @@ def __init__( parameters: dict[str, float], quadrature_rule: QuadratureRule, q_fields_local: dict[str, df.fem.Function], - q_field_nonlocal_rate: df.fem.Function, + q_field_nonlocal: df.fem.Function, Q_local_damage: str | None = None, ): q_fields = {Q_local: q_fields_local[Q_local]} @@ -387,9 +387,8 @@ def __init__( if Q_nonlocal_rate[-4:] != "rate": raise ValueError("Q_nonlocal_rate must be some rate for CDM, you provided " + Q_nonlocal_rate) - q_fields[Q_nonlocal_rate] = q_field_nonlocal_rate - Q_nonlocal = Q_nonlocal_rate[:-5] + q_fields[Q_nonlocal] = q_field_nonlocal fields = { Q_nonlocal_rate: df.fem.Function(function_space, name=Q_nonlocal_rate), @@ -419,7 +418,7 @@ def __init__( f_form = df.fem.form(f_ufl) - rate_evaluator = QuadratureEvaluator(fields[Q_nonlocal_rate], function_space.mesh, quadrature_rule) + strain_evaluator = QuadratureEvaluator(fields[Q_nonlocal], function_space.mesh, quadrature_rule) super().__init__( Q_local=Q_local, Q_nonlocal=Q_nonlocal, @@ -431,7 +430,7 @@ def __init__( q_fields=q_fields, fields=fields, form=f_form, - rate_evaluator=rate_evaluator, + strain_evaluator=strain_evaluator, ) def step(self, h: float) -> None: @@ -454,8 +453,8 @@ def step(self, h: float) -> None: self.fields[self.Q_nonlocal].vector.array[:] += h * self.fields[self.Q_nonlocal_rate].vector.array self.fields[self.Q_nonlocal].x.scatter_forward() - self.rate_evaluator(self.q_fields[self.Q_nonlocal_rate]) - self.q_fields[self.Q_nonlocal_rate].x.scatter_forward() + self.strain_evaluator(self.q_fields[self.Q_nonlocal]) + self.q_fields[self.Q_nonlocal].x.scatter_forward() self.t += h diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index a02efdb..bbf21a3 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -51,7 +51,12 @@ impl ConstitutiveModel for GradientJH23D { d_eps_vol *= -1.0; let sigma_0 = input.get_vector::<{ Q::MandelStress.size() }>(Q::MandelStress, ip); - let del_lambda_nonlocal = del_t * input.get_scalar(Q::EqNonlocalPlasticStrainRate, ip).max(0.0); + //let del_lambda_nonlocal = del_t * input.get_scalar(Q::EqNonlocalPlasticStrainRate, ip).max(0.0); + let del_lambda_nonlocal = (input.get_scalar(Q::HistoryMaximum, ip) - input.get_scalar(Q::EqNonlocalPlasticStrain, ip)).max(0.0); + if del_lambda_nonlocal > 0.0 { + output.set_scalar(Q::HistoryMaximum, ip, input.get_scalar(Q::EqNonlocalPlasticStrain, ip)); + } + let mut del_lambda = 0.0; @@ -74,9 +79,7 @@ impl ConstitutiveModel for GradientJH23D { let damage_0 = input.get_scalar(Q::Damage, ip); let mut damage_1 = 0.0; if self.parameters.E_F > 0.0 { - let lambda_old = - self.parameters.E_F * (1. - damage_0).ln(); - let lambda_new = lambda_old + del_lambda_nonlocal; - damage_1 = 1. - (-lambda_new / self.parameters.E_F).exp(); + damage_1 = 1. - (-input.get_scalar(Q::EqNonlocalPlasticStrain, ip) / self.parameters.E_F).exp(); } else { damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); } @@ -229,7 +232,8 @@ impl ConstitutiveModel for GradientJH23D { fn define_input(&self) -> HashMap { HashMap::from([ (Q::VelocityGradient, QDim::SquareTensor(3)), - (Q::EqNonlocalPlasticStrainRate, QDim::Scalar), + //(Q::EqNonlocalPlasticStrainRate, QDim::Scalar), + (Q::EqNonlocalPlasticStrain, QDim::Scalar), (Q::CellDiameter, QDim::Scalar), ]) } @@ -244,6 +248,7 @@ impl ConstitutiveModel for GradientJH23D { (Q::Density, QDim::Scalar), (Q::EqPlasticStrain, QDim::Scalar), (Q::BulkViscosity, QDim::Scalar), + (Q::HistoryMaximum, QDim::Scalar), ]) } diff --git a/src/interfaces.rs b/src/interfaces.rs index c742c1c..7abf451 100644 --- a/src/interfaces.rs +++ b/src/interfaces.rs @@ -63,6 +63,8 @@ pub enum Q { BulkViscosity, #[strum(serialize = "CellDiameter", serialize = "cell_diameter")] CellDiameter, + #[strum(serialize = "HistoryMaximum", serialize = "history_maximum")] + HistoryMaximum, #[strum(serialize = "_LAST", serialize = "_last")] _LAST, } @@ -234,6 +236,7 @@ impl Q { Q::InternalPlasticEnergyRate => QDim::Scalar, Q::BulkViscosity => QDim::Scalar, Q::CellDiameter => QDim::Scalar, + Q::HistoryMaximum => QDim::Scalar, Q::_LAST => QDim::Scalar, } } From 2ba814d9ba70b2aebe6d06b7edadc104014d08c8 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 14 Feb 2024 15:23:17 +0100 Subject: [PATCH 06/24] fix some errors --- python/comfe/cdm.py | 2 +- src/gradient_jh2.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/python/comfe/cdm.py b/python/comfe/cdm.py index 5a1c952..daf3519 100644 --- a/python/comfe/cdm.py +++ b/python/comfe/cdm.py @@ -678,7 +678,7 @@ def __init__( nonlocal_parameters, quadrature_rule, mechanics_solver.q_fields, - mechanics_solver.model.input[Q_nonlocal_rate], + mechanics_solver.model.input[Q_nonlocal], Q_local_damage=Q_local_damage, ) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index bbf21a3..d2005a1 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -52,7 +52,7 @@ impl ConstitutiveModel for GradientJH23D { let sigma_0 = input.get_vector::<{ Q::MandelStress.size() }>(Q::MandelStress, ip); //let del_lambda_nonlocal = del_t * input.get_scalar(Q::EqNonlocalPlasticStrainRate, ip).max(0.0); - let del_lambda_nonlocal = (input.get_scalar(Q::HistoryMaximum, ip) - input.get_scalar(Q::EqNonlocalPlasticStrain, ip)).max(0.0); + let del_lambda_nonlocal = (input.get_scalar(Q::EqNonlocalPlasticStrain, ip)- input.get_scalar(Q::HistoryMaximum, ip)).max(0.0); if del_lambda_nonlocal > 0.0 { output.set_scalar(Q::HistoryMaximum, ip, input.get_scalar(Q::EqNonlocalPlasticStrain, ip)); } From 0247a60194dd81892a8da715b9d90abb4507851a Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 14 Feb 2024 17:04:30 +0100 Subject: [PATCH 07/24] fix another error --- src/gradient_jh2.rs | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index d2005a1..ee28460 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -52,10 +52,12 @@ impl ConstitutiveModel for GradientJH23D { let sigma_0 = input.get_vector::<{ Q::MandelStress.size() }>(Q::MandelStress, ip); //let del_lambda_nonlocal = del_t * input.get_scalar(Q::EqNonlocalPlasticStrainRate, ip).max(0.0); - let del_lambda_nonlocal = (input.get_scalar(Q::EqNonlocalPlasticStrain, ip)- input.get_scalar(Q::HistoryMaximum, ip)).max(0.0); - if del_lambda_nonlocal > 0.0 { - output.set_scalar(Q::HistoryMaximum, ip, input.get_scalar(Q::EqNonlocalPlasticStrain, ip)); - } + let nonlocal_plastic_strain = input.get_scalar(Q::EqNonlocalPlasticStrain, ip); + let history_maximum_0 = input.get_scalar(Q::HistoryMaximum, ip); + //let del_lambda_nonlocal = input.get_scalar(Q::EqNonlocalPlasticStrain, ip)- input.get_scalar(Q::HistoryMaximum, ip); + let history_maximum_1 = history_maximum_0.max(nonlocal_plastic_strain); + let del_lambda_nonlocal = history_maximum_1 - history_maximum_0; + output.set_scalar(Q::HistoryMaximum, ip, history_maximum_1); let mut del_lambda = 0.0; @@ -79,7 +81,7 @@ impl ConstitutiveModel for GradientJH23D { let damage_0 = input.get_scalar(Q::Damage, ip); let mut damage_1 = 0.0; if self.parameters.E_F > 0.0 { - damage_1 = 1. - (-input.get_scalar(Q::EqNonlocalPlasticStrain, ip) / self.parameters.E_F).exp(); + damage_1 = 1. - (-history_maximum_1 / self.parameters.E_F).exp(); } else { damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); } From 234bc7689efc48117b08bfed13817ae96d3dcf61 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Mon, 19 Feb 2024 11:08:55 +0100 Subject: [PATCH 08/24] WIP: nonlocal balance in init conf --- python/comfe/cdm.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/python/comfe/cdm.py b/python/comfe/cdm.py index daf3519..268221e 100644 --- a/python/comfe/cdm.py +++ b/python/comfe/cdm.py @@ -381,6 +381,7 @@ def __init__( q_fields_local: dict[str, df.fem.Function], q_field_nonlocal: df.fem.Function, Q_local_damage: str | None = None, + density_0: float | None = None, ): q_fields = {Q_local: q_fields_local[Q_local]} @@ -407,12 +408,21 @@ def __init__( else: g = 1.0 + if density_0 is not None: + frac_det_F = q_fields["density"] * (1.0 / density_0) + else: + frac_det_F = 1.0 + f_int_ufl = ( - g * parameters["l"] ** 2 * ufl.inner(ufl.grad(fields[Q_nonlocal]), ufl.grad(test_function)) - + fields[Q_nonlocal] * test_function - ) * quadrature_rule.dx + frac_det_F + * ( + g * parameters["l"] ** 2 * ufl.inner(ufl.grad(fields[Q_nonlocal]), ufl.grad(test_function)) + + fields[Q_nonlocal] * test_function + ) + * quadrature_rule.dx + ) - f_ext_ufl = q_fields[Q_local] * test_function * quadrature_rule.dx + f_ext_ufl = frac_det_F * q_fields[Q_local] * test_function * quadrature_rule.dx f_ufl = -f_int_ufl + f_ext_ufl @@ -680,6 +690,7 @@ def __init__( mechanics_solver.q_fields, mechanics_solver.model.input[Q_nonlocal], Q_local_damage=Q_local_damage, + density_0=parameters["rho"], ) # add all fields from the solver to this class for easier postprocessing From 803e47307c0a680ba083f0dde0d253c7505c1b45 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Mon, 19 Feb 2024 11:36:09 +0100 Subject: [PATCH 09/24] fix key error --- python/comfe/cdm.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/comfe/cdm.py b/python/comfe/cdm.py index 268221e..d0619aa 100644 --- a/python/comfe/cdm.py +++ b/python/comfe/cdm.py @@ -409,6 +409,7 @@ def __init__( g = 1.0 if density_0 is not None: + q_fields["density"] = q_fields_local["density"] frac_det_F = q_fields["density"] * (1.0 / density_0) else: frac_det_F = 1.0 From 3f2abe36aadd0781705469668db2e600d3c661bc Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 21 Feb 2024 12:25:05 +0100 Subject: [PATCH 10/24] initial_config in implicit --- python/comfe/cdm.py | 46 +++++++++++++++++++++++++++++---------------- 1 file changed, 30 insertions(+), 16 deletions(-) diff --git a/python/comfe/cdm.py b/python/comfe/cdm.py index d0619aa..4347db5 100644 --- a/python/comfe/cdm.py +++ b/python/comfe/cdm.py @@ -474,7 +474,7 @@ class ImplicitNonlocal(NonlocalInterface): # M: df.fem.Function # fields: dict[str, df.fem.Function] # q_fields: dict[str, df.fem.Function | ufl.core.expr.Expr] - rate_evaluator: QuadratureEvaluator + strain_evaluator: QuadratureEvaluator parameters: dict[str, float] # form: df.fem.FormMetaClass problem: df.fem.petsc.LinearProblem @@ -493,20 +493,22 @@ def __init__( parameters: dict[str, float], quadrature_rule: QuadratureRule, q_fields_local: dict[str, df.fem.Function], - q_field_nonlocal_rate: df.fem.Function, + q_field_nonlocal: df.fem.Function, Q_local_damage: str | None = None, + density_0: float | None = None, ): q_fields = {Q_local: q_fields_local[Q_local]} if Q_nonlocal_rate[-4:] != "rate": raise ValueError("Q_nonlocal_rate must be some rate for CDM, you provided " + Q_nonlocal_rate) - q_fields[Q_nonlocal_rate] = q_field_nonlocal_rate + # q_fields[Q_nonlocal_rate] = q_field_nonlocal_rate Q_nonlocal = Q_nonlocal_rate[:-5] + q_fields[Q_nonlocal] = q_field_nonlocal fields = { - Q_nonlocal_rate: df.fem.Function(function_space, name=Q_nonlocal_rate), + # Q_nonlocal_rate: df.fem.Function(function_space, name=Q_nonlocal_rate), Q_nonlocal: df.fem.Function(function_space, name=Q_nonlocal), "nonlocal_force": df.fem.Function(function_space, name="nonlocal_force"), } @@ -520,15 +522,25 @@ def __init__( else: g = 1.0 + if density_0 is not None: + q_fields["density"] = q_fields_local["density"] + frac_det_F = q_fields["density"] * (1.0 / density_0) + else: + frac_det_F = 1.0 + test_function = ufl.TestFunction(function_space) trial_function = ufl.TrialFunction(function_space) - b_form = ufl.inner(test_function, q_fields[Q_local]) * quadrature_rule.dx + b_form = frac_det_F * ufl.inner(test_function, q_fields[Q_local]) * quadrature_rule.dx A_form = ( - ufl.inner(test_function, trial_function) - + g * parameters["l"] ** 2.0 * ufl.inner(ufl.grad(test_function), ufl.grad(trial_function)) - ) * quadrature_rule.dx + frac_det_F + * ( + ufl.inner(test_function, trial_function) + + g * parameters["l"] ** 2.0 * ufl.inner(ufl.grad(test_function), ufl.grad(trial_function)) + ) + * quadrature_rule.dx + ) problem = df.fem.petsc.LinearProblem(A_form, b_form, u=fields[Q_nonlocal]) @@ -543,7 +555,7 @@ def __init__( # f_form = df.fem.form(f_ufl) - rate_evaluator = QuadratureEvaluator(fields[Q_nonlocal_rate], function_space.mesh, quadrature_rule) + strain_evaluator = QuadratureEvaluator(fields[Q_nonlocal], function_space.mesh, quadrature_rule) super().__init__( Q_local=Q_local, Q_nonlocal=Q_nonlocal, @@ -555,17 +567,17 @@ def __init__( q_fields=q_fields, fields=fields, # form=f_form, - rate_evaluator=rate_evaluator, + strain_evaluator=strain_evaluator, problem=problem, ) def step(self, h: float) -> None: - old_nonlocal = self.fields[self.Q_nonlocal].vector.array.copy() + # old_nonlocal = self.fields[self.Q_nonlocal].vector.array.copy() self.problem.solve() - self.q_fields[self.Q_nonlocal_rate].vector.array[:] = ( - self.fields[self.Q_nonlocal].vector.array - old_nonlocal - ) / h - self.rate_evaluator(self.q_fields[self.Q_nonlocal_rate]) + # self.q_fields[self.Q_nonlocal_rate].vector.array[:] = ( + # self.fields[self.Q_nonlocal].vector.array - old_nonlocal + # ) / h + self.strain_evaluator(self.q_fields[self.Q_nonlocal]) self.t += h @@ -657,6 +669,7 @@ def __init__( mass_mechanics: df.fem.Function | None = None, mass_nonlocal: df.fem.Function | None = None, calculate_bulk_viscosity: bool = False, + nonlocal_initial_config: bool = True, ) -> None: mass_mechanics = ( mass_mechanics @@ -680,6 +693,7 @@ def __init__( calculate_bulk_viscosity=calculate_bulk_viscosity, ) + density_0 = parameters["rho"] if nonlocal_initial_config else None nonlocal_solver = nonlocal_solver( Q_local, Q_nonlocal_rate, @@ -691,7 +705,7 @@ def __init__( mechanics_solver.q_fields, mechanics_solver.model.input[Q_nonlocal], Q_local_damage=Q_local_damage, - density_0=parameters["rho"], + density_0=density_0, ) # add all fields from the solver to this class for easier postprocessing From 73596275dbade3d515a20d5ab784ec38217f99ad Mon Sep 17 00:00:00 2001 From: srosenbu Date: Mon, 26 Feb 2024 13:13:12 +0100 Subject: [PATCH 11/24] Fix error in initial connfig --- python/comfe/cdm.py | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/python/comfe/cdm.py b/python/comfe/cdm.py index 4347db5..3d82c61 100644 --- a/python/comfe/cdm.py +++ b/python/comfe/cdm.py @@ -381,7 +381,7 @@ def __init__( q_fields_local: dict[str, df.fem.Function], q_field_nonlocal: df.fem.Function, Q_local_damage: str | None = None, - density_0: float | None = None, + displacements: df.fem.Function | None = None, ): q_fields = {Q_local: q_fields_local[Q_local]} @@ -408,22 +408,17 @@ def __init__( else: g = 1.0 - if density_0 is not None: - q_fields["density"] = q_fields_local["density"] - frac_det_F = q_fields["density"] * (1.0 / density_0) + if displacements is not None: + fields["u"] = displacements else: - frac_det_F = 1.0 + fields["u"] = None f_int_ufl = ( - frac_det_F - * ( - g * parameters["l"] ** 2 * ufl.inner(ufl.grad(fields[Q_nonlocal]), ufl.grad(test_function)) - + fields[Q_nonlocal] * test_function - ) - * quadrature_rule.dx - ) + g * parameters["l"] ** 2 * ufl.inner(ufl.grad(fields[Q_nonlocal]), ufl.grad(test_function)) + + fields[Q_nonlocal] * test_function + ) * quadrature_rule.dx - f_ext_ufl = frac_det_F * q_fields[Q_local] * test_function * quadrature_rule.dx + f_ext_ufl = q_fields[Q_local] * test_function * quadrature_rule.dx f_ufl = -f_int_ufl + f_ext_ufl @@ -448,6 +443,9 @@ def step(self, h: float) -> None: with self.fields["nonlocal_force"].vector.localForm() as f_local: f_local.set(0.0) + if self.fields["u"] is not None: + self.fields["u"].function_space.mesh.geometry.x[:] -= self.fields["u"].x.array + df.fem.petsc.assemble_vector(self.fields["nonlocal_force"].vector, self.form) self.fields["nonlocal_force"].x.scatter_reverse(ScatterMode.add) @@ -467,6 +465,9 @@ def step(self, h: float) -> None: self.strain_evaluator(self.q_fields[self.Q_nonlocal]) self.q_fields[self.Q_nonlocal].x.scatter_forward() + if self.fields["u"] is not None: + self.fields["u"].function_space.mesh.geometry.x[:] += self.fields["u"].x.array + self.t += h @@ -693,7 +694,7 @@ def __init__( calculate_bulk_viscosity=calculate_bulk_viscosity, ) - density_0 = parameters["rho"] if nonlocal_initial_config else None + displacements = mechanics_solver.fields["u"] if nonlocal_initial_config else None nonlocal_solver = nonlocal_solver( Q_local, Q_nonlocal_rate, @@ -705,7 +706,7 @@ def __init__( mechanics_solver.q_fields, mechanics_solver.model.input[Q_nonlocal], Q_local_damage=Q_local_damage, - density_0=density_0, + displacements=displacements, ) # add all fields from the solver to this class for easier postprocessing From bdfc7ee58dd350f7231c2fb6c4b44551de2c2fdd Mon Sep 17 00:00:00 2001 From: srosenbu Date: Mon, 26 Feb 2024 13:19:39 +0100 Subject: [PATCH 12/24] error in mesh update --- python/comfe/cdm.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/python/comfe/cdm.py b/python/comfe/cdm.py index 3d82c61..bd7adcf 100644 --- a/python/comfe/cdm.py +++ b/python/comfe/cdm.py @@ -444,7 +444,8 @@ def step(self, h: float) -> None: f_local.set(0.0) if self.fields["u"] is not None: - self.fields["u"].function_space.mesh.geometry.x[:] -= self.fields["u"].x.array + set_mesh_coordinates(self.fields["u"].function_space.mesh, -self.fields["u"].x.array, mode="add") + # self.fields["u"].function_space.mesh.geometry.x[:] -= self.fields["u"].x.array df.fem.petsc.assemble_vector(self.fields["nonlocal_force"].vector, self.form) self.fields["nonlocal_force"].x.scatter_reverse(ScatterMode.add) @@ -466,7 +467,8 @@ def step(self, h: float) -> None: self.q_fields[self.Q_nonlocal].x.scatter_forward() if self.fields["u"] is not None: - self.fields["u"].function_space.mesh.geometry.x[:] += self.fields["u"].x.array + set_mesh_coordinates(self.fields["u"].function_space.mesh, self.fields["u"].x.array, mode="add") + # self.fields["u"].function_space.mesh.geometry.x[:] += self.fields["u"].x.array self.t += h From 108ddae31e760dfc4b89d2100cbcd2f3b9e1534a Mon Sep 17 00:00:00 2001 From: srosenbu Date: Mon, 26 Feb 2024 17:33:31 +0100 Subject: [PATCH 13/24] fix condition for init config --- python/comfe/cdm.py | 28 +++++++++++----------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/python/comfe/cdm.py b/python/comfe/cdm.py index bd7adcf..aaa6ba6 100644 --- a/python/comfe/cdm.py +++ b/python/comfe/cdm.py @@ -410,8 +410,6 @@ def __init__( if displacements is not None: fields["u"] = displacements - else: - fields["u"] = None f_int_ufl = ( g * parameters["l"] ** 2 * ufl.inner(ufl.grad(fields[Q_nonlocal]), ufl.grad(test_function)) @@ -443,7 +441,7 @@ def step(self, h: float) -> None: with self.fields["nonlocal_force"].vector.localForm() as f_local: f_local.set(0.0) - if self.fields["u"] is not None: + if "u" in self.fields: set_mesh_coordinates(self.fields["u"].function_space.mesh, -self.fields["u"].x.array, mode="add") # self.fields["u"].function_space.mesh.geometry.x[:] -= self.fields["u"].x.array @@ -466,7 +464,7 @@ def step(self, h: float) -> None: self.strain_evaluator(self.q_fields[self.Q_nonlocal]) self.q_fields[self.Q_nonlocal].x.scatter_forward() - if self.fields["u"] is not None: + if "u" in self.fields: set_mesh_coordinates(self.fields["u"].function_space.mesh, self.fields["u"].x.array, mode="add") # self.fields["u"].function_space.mesh.geometry.x[:] += self.fields["u"].x.array @@ -498,7 +496,7 @@ def __init__( q_fields_local: dict[str, df.fem.Function], q_field_nonlocal: df.fem.Function, Q_local_damage: str | None = None, - density_0: float | None = None, + displacements: df.fem.Function | None = None, ): q_fields = {Q_local: q_fields_local[Q_local]} @@ -525,25 +523,20 @@ def __init__( else: g = 1.0 - if density_0 is not None: - q_fields["density"] = q_fields_local["density"] - frac_det_F = q_fields["density"] * (1.0 / density_0) + if displacements is not None: + fields["u"] = displacements else: - frac_det_F = 1.0 + fields["u"] = None test_function = ufl.TestFunction(function_space) trial_function = ufl.TrialFunction(function_space) - b_form = frac_det_F * ufl.inner(test_function, q_fields[Q_local]) * quadrature_rule.dx + b_form = ufl.inner(test_function, q_fields[Q_local]) * quadrature_rule.dx A_form = ( - frac_det_F - * ( - ufl.inner(test_function, trial_function) - + g * parameters["l"] ** 2.0 * ufl.inner(ufl.grad(test_function), ufl.grad(trial_function)) - ) - * quadrature_rule.dx - ) + ufl.inner(test_function, trial_function) + + g * parameters["l"] ** 2.0 * ufl.inner(ufl.grad(test_function), ufl.grad(trial_function)) + ) * quadrature_rule.dx problem = df.fem.petsc.LinearProblem(A_form, b_form, u=fields[Q_nonlocal]) @@ -576,6 +569,7 @@ def __init__( def step(self, h: float) -> None: # old_nonlocal = self.fields[self.Q_nonlocal].vector.array.copy() + self.problem.solve() # self.q_fields[self.Q_nonlocal_rate].vector.array[:] = ( # self.fields[self.Q_nonlocal].vector.array - old_nonlocal From adf75c9de05160c0a2044d71453b9f4166d7fb06 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Tue, 12 Mar 2024 12:35:39 +0100 Subject: [PATCH 14/24] Fix --- python/comfe/cdm.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/python/comfe/cdm.py b/python/comfe/cdm.py index aaa6ba6..03ff37b 100644 --- a/python/comfe/cdm.py +++ b/python/comfe/cdm.py @@ -525,8 +525,6 @@ def __init__( if displacements is not None: fields["u"] = displacements - else: - fields["u"] = None test_function = ufl.TestFunction(function_space) trial_function = ufl.TrialFunction(function_space) @@ -569,12 +567,16 @@ def __init__( def step(self, h: float) -> None: # old_nonlocal = self.fields[self.Q_nonlocal].vector.array.copy() + if "u" in self.fields: + set_mesh_coordinates(self.fields["u"].function_space.mesh, -self.fields["u"].x.array, mode="add") self.problem.solve() # self.q_fields[self.Q_nonlocal_rate].vector.array[:] = ( # self.fields[self.Q_nonlocal].vector.array - old_nonlocal # ) / h self.strain_evaluator(self.q_fields[self.Q_nonlocal]) + if "u" in self.fields: + set_mesh_coordinates(self.fields["u"].function_space.mesh, self.fields["u"].x.array, mode="add") self.t += h @@ -675,7 +677,7 @@ def __init__( ) mass_nonlocal = ( mass_nonlocal - if mass_nonlocal is not None or nonlocal_parameters["solver"] == "implicit" + if mass_nonlocal is not None or nonlocal_solver == ImplicitNonlocal else diagonal_mass(nonlocal_space, nonlocal_parameters["zeta"], invert=True) ) mechanics_solver = mechanics_solver( From 25bf38a78a61b9479345bc576439f88b194c4602 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Tue, 12 Mar 2024 12:36:57 +0100 Subject: [PATCH 15/24] local sound speed --- src/gradient_jh2.rs | 13 ++++++++++--- src/jh2.rs | 29 +++++++++++++++++------------ 2 files changed, 27 insertions(+), 15 deletions(-) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index ee28460..d950378 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -36,6 +36,8 @@ impl ConstitutiveModel for GradientJH23D { EFMIN: *parameters.get("EFMIN").unwrap(), DMAX: *parameters.get("DMAX").unwrap_or(&1.0), E_F: *parameters.get("E_F").unwrap_or(&0.0), + LOCAL_SOUND_SPEED: *parameters.get("LOCAL_SOUND_SPEED").unwrap_or(&0.0), + //M_OVERNONLOCAL: *parameters.get("M_OVERNONLOCAL").unwrap_or(&0.0), //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, }) @@ -125,19 +127,22 @@ impl ConstitutiveModel for GradientJH23D { output.set_scalar(Q::Density, ip, density_1); let mu = density_1 / self.parameters.RHO - 1.; - + let mut local_bulk_modulus = self.parameters.K1; let p_1 = { if mu > 0.0 { + local_bulk_modulus += + 2.0 * self.parameters.K2 * mu + 3.0 * self.parameters.K3 * mu.powi(2); self.parameters.K1 * mu + self.parameters.K2 * mu.powi(2) + self.parameters.K3 * mu.powi(3) + input.get_scalar(Q::BulkingPressure, ip) } else { - let p_trial = self.parameters.K1*mu; + let p_trial = self.parameters.K1 * mu; let p_damaged = -self.parameters.T * (1. - damage_1); if p_trial > p_damaged { p_trial } else { + local_bulk_modulus = 0.0; p_damaged } } @@ -161,7 +166,9 @@ impl ConstitutiveModel for GradientJH23D { // Calculate bulk viscosity let l = input.get_scalar(Q::CellDiameter, ip); let tr_deps = d_eps_vol * 3.; - let c = (self.parameters.K1 / self.parameters.RHO).sqrt(); + let bulk_modulus = local_bulk_modulus * self.parameters.LOCAL_SOUND_SPEED + + (1.0 - self.parameters.LOCAL_SOUND_SPEED) * self.parameters.K1; + let c = (bulk_modulus / self.parameters.RHO).sqrt(); let q_1 = { if tr_deps < 0.0 { l * density_1 * (1.5 * l * tr_deps.powi(2) + c * 0.06 * tr_deps.abs()) diff --git a/src/jh2.rs b/src/jh2.rs index 22dfd8e..b90370b 100644 --- a/src/jh2.rs +++ b/src/jh2.rs @@ -26,6 +26,7 @@ pub struct JH2ConstParameters { pub EFMIN: f64, pub DMAX: f64, pub E_F: f64, + pub LOCAL_SOUND_SPEED: f64, //pub REDUCE_T: f64, } #[derive(Debug)] @@ -57,6 +58,7 @@ impl ConstitutiveModel for JH23D { EFMIN: *parameters.get("EFMIN").unwrap(), DMAX: *parameters.get("DMAX").unwrap_or(&1.0), E_F: *parameters.get("E_F").unwrap_or(&0.0), + LOCAL_SOUND_SPEED: *parameters.get("LOCAL_SOUND_SPEED").unwrap_or(&0.0), //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, }) @@ -89,7 +91,7 @@ impl ConstitutiveModel for JH23D { let p_s = p_0 / self.parameters.PHEL; let t_s = self.parameters.T / self.parameters.PHEL; let mut rate_factor = 1.; - + //let t_s_factor = (1.-damage_0).powf(self.parameters.REDUCE_T); let fracture_surface = (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) @@ -114,7 +116,7 @@ impl ConstitutiveModel for JH23D { alpha = yield_surface / s_tr_eq; if self.parameters.E_F > 0.0 { - let lambda_old = - self.parameters.E_F * (1. - damage_0).ln(); + let lambda_old = -self.parameters.E_F * (1. - damage_0).ln(); let lambda_new = lambda_old + del_lambda; damage_1 = 1. - (-lambda_new / self.parameters.E_F).exp(); } else { @@ -143,9 +145,11 @@ impl ConstitutiveModel for JH23D { output.set_scalar(Q::Density, ip, density_1); let mu = density_1 / self.parameters.RHO - 1.; - + let mut local_bulk_modulus = self.parameters.K1; let p_1 = { if mu > 0.0 { + local_bulk_modulus += + 2.0 * self.parameters.K2 * mu + 3.0 * self.parameters.K3 * mu.powi(2); self.parameters.K1 * mu + self.parameters.K2 * mu.powi(2) + self.parameters.K3 * mu.powi(3) @@ -156,6 +160,7 @@ impl ConstitutiveModel for JH23D { if p_trial > p_damaged { p_trial } else { + local_bulk_modulus = 0.0; p_damaged } } @@ -179,7 +184,9 @@ impl ConstitutiveModel for JH23D { // Calculate bulk viscosity let l = input.get_scalar(Q::CellDiameter, ip); let tr_deps = d_eps_vol * 3.; - let c = (self.parameters.K1 / self.parameters.RHO).sqrt(); + let bulk_modulus = local_bulk_modulus * self.parameters.LOCAL_SOUND_SPEED + + (1.0 - self.parameters.LOCAL_SOUND_SPEED) * self.parameters.K1; + let c = (bulk_modulus / self.parameters.RHO).sqrt(); let q_1 = { if tr_deps < 0.0 { l * density_1 * (1.5 * l * tr_deps.powi(2) + c * 0.06 * tr_deps.abs()) @@ -213,15 +220,14 @@ impl ConstitutiveModel for JH23D { } let elastic_rate = - - (1. - alpha) / (2. * self.parameters.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; - + -(1. - alpha) / (2. * self.parameters.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; + let density_mid = 0.5 * (density_0 + density_1); if output.is_some(Q::InternalPlasticEnergy) && input.is_some(Q::InternalPlasticEnergy) { let s_mid = 0.5 * (s_0 + s_1); let deviatoric_rate = d_eps_dev - elastic_rate; let e_0 = input.get_scalar(Q::InternalPlasticEnergy, ip); - let e_1 = e_0 - + del_t / density_mid * (s_mid.dot(&deviatoric_rate)); + let e_1 = e_0 + del_t / density_mid * (s_mid.dot(&deviatoric_rate)); output.set_scalar(Q::InternalPlasticEnergy, ip, e_1); } if output.is_some(Q::InternalElasticEnergy) && input.is_some(Q::InternalElasticEnergy) { @@ -229,9 +235,8 @@ impl ConstitutiveModel for JH23D { let p_mid = -0.5 * (p_0 + p_1); let deviatoric_rate = elastic_rate; let e_0 = input.get_scalar(Q::InternalElasticEnergy, ip); - let e_1 = e_0 - + del_t / density_mid - * (s_mid.dot(&deviatoric_rate) + 3. * d_eps_vol * p_mid); + let e_1 = + e_0 + del_t / density_mid * (s_mid.dot(&deviatoric_rate) + 3. * d_eps_vol * p_mid); output.set_scalar(Q::InternalElasticEnergy, ip, e_1); } if output.is_some(Q::InternalEnergy) && input.is_some(Q::InternalEnergy) { @@ -242,7 +247,7 @@ impl ConstitutiveModel for JH23D { } if output.is_some(Q::InternalHeatingEnergy) && input.is_some(Q::InternalHeatingEnergy) { let e_0 = input.get_scalar(Q::InternalHeatingEnergy, ip); - let q_mid = - 0.5 * (q_1 + input.get_scalar(Q::BulkViscosity, ip)); + let q_mid = -0.5 * (q_1 + input.get_scalar(Q::BulkViscosity, ip)); let e_1 = e_0 + del_t / density_mid * 3. * q_mid * d_eps_vol; output.set_scalar(Q::InternalHeatingEnergy, ip, e_1); } From ecfc3da43327924e769616b78c42f0a34125c78d Mon Sep 17 00:00:00 2001 From: srosenbu Date: Tue, 12 Mar 2024 12:37:29 +0100 Subject: [PATCH 16/24] stuff --- test/test_jh2_2d.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/test/test_jh2_2d.py b/test/test_jh2_2d.py index adc7659..c72ee92 100644 --- a/test/test_jh2_2d.py +++ b/test/test_jh2_2d.py @@ -503,11 +503,6 @@ def test_single_element_2d(test_case: dict, plot: str | None = None) -> None: p_ = np.array(p_).reshape((-1, 1)) s_eq_ = np.array(s_eq_).reshape((-1, 1)) - points = np.hstack((p_, s_eq_)) - tree = KDTree(points) - distances = tree.query(test_case["points"]) - assert np.mean(distances[0] / np.max(np.abs(test_case["points"][:, 1]))) < 0.05 - if plot is not None: import matplotlib import matplotlib.pyplot as plt @@ -524,7 +519,13 @@ def test_single_element_2d(test_case: dict, plot: str | None = None) -> None: plt.savefig(f"{plot}.png") plt.clf() + points = np.hstack((p_, s_eq_)) + tree = KDTree(points) + distances = tree.query(test_case["points"]) + assert np.mean(distances[0] / np.max(np.abs(test_case["points"][:, 1]))) < 0.05 + if __name__ == "__main__": test_single_element_2d(case_3a, plot="jh2_case_3a") test_single_element_2d(case_1, plot="jh2_case_1") + test_single_element_2d(case_2, plot="jh2_case_2") From ce8ea4df9a7ba66cd3cb23067cfbe698793c7872 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 13 Mar 2024 10:19:26 +0100 Subject: [PATCH 17/24] Overnonlocal model --- src/gradient_jh2.rs | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index d950378..a4882c7 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -7,9 +7,15 @@ use crate::jh2::JH2ConstParameters; use std::collections::HashMap; +#[derive(Debug)] +pub struct NonlocalParameters{ + m_overnonlocal: f64, +} + #[derive(Debug)] pub struct GradientJH23D { parameters: JH2ConstParameters, + nonlocal_parameters: NonlocalParameters, } impl ConstitutiveModel for GradientJH23D { @@ -40,6 +46,9 @@ impl ConstitutiveModel for GradientJH23D { //M_OVERNONLOCAL: *parameters.get("M_OVERNONLOCAL").unwrap_or(&0.0), //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, + nonlocal_parameters: NonlocalParameters { + m_overnonlocal: *parameters.get("M_OVERNONLOCAL").unwrap_or(&1.0), + }, }) } fn evaluate_ip(&self, ip: usize, del_t: f64, input: &QValueInput, output: &mut QValueOutput) { @@ -83,7 +92,9 @@ impl ConstitutiveModel for GradientJH23D { let damage_0 = input.get_scalar(Q::Damage, ip); let mut damage_1 = 0.0; if self.parameters.E_F > 0.0 { - damage_1 = 1. - (-history_maximum_1 / self.parameters.E_F).exp(); + let m = self.nonlocal_parameters.m_overnonlocal; + let kappa = m * history_maximum_1 + (1.0 - m) * input.get_scalar(Q::EqPlasticStrain, ip); + damage_1 = 1. - (-kappa / self.parameters.E_F).exp(); } else { damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); } From 474a5bdd15b54c5e8d12027ad73e6f0173032669 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 13 Mar 2024 13:02:03 +0100 Subject: [PATCH 18/24] Fix negative damage --- src/gradient_jh2.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index a4882c7..a650326 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -93,7 +93,7 @@ impl ConstitutiveModel for GradientJH23D { let mut damage_1 = 0.0; if self.parameters.E_F > 0.0 { let m = self.nonlocal_parameters.m_overnonlocal; - let kappa = m * history_maximum_1 + (1.0 - m) * input.get_scalar(Q::EqPlasticStrain, ip); + let kappa = (m * history_maximum_1 + (1.0 - m) * input.get_scalar(Q::EqPlasticStrain, ip)).max(0.0); damage_1 = 1. - (-kappa / self.parameters.E_F).exp(); } else { damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); From 17a122667a50ce730dded45f4e715b1df775e60f Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 3 Apr 2024 11:21:53 +0200 Subject: [PATCH 19/24] e_0 in damage --- src/gradient_jh2.rs | 4 +++- src/jh2.rs | 6 +++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index a650326..40d9772 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -42,7 +42,9 @@ impl ConstitutiveModel for GradientJH23D { EFMIN: *parameters.get("EFMIN").unwrap(), DMAX: *parameters.get("DMAX").unwrap_or(&1.0), E_F: *parameters.get("E_F").unwrap_or(&0.0), + E_0: *parameters.get("E_0").unwrap_or(&0.0), LOCAL_SOUND_SPEED: *parameters.get("LOCAL_SOUND_SPEED").unwrap_or(&0.0), + HARDENING_MODULUS: *parameters.get("HARDENING_MODULUS").unwrap_or(&0.0), //M_OVERNONLOCAL: *parameters.get("M_OVERNONLOCAL").unwrap_or(&0.0), //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, @@ -94,7 +96,7 @@ impl ConstitutiveModel for GradientJH23D { if self.parameters.E_F > 0.0 { let m = self.nonlocal_parameters.m_overnonlocal; let kappa = (m * history_maximum_1 + (1.0 - m) * input.get_scalar(Q::EqPlasticStrain, ip)).max(0.0); - damage_1 = 1. - (-kappa / self.parameters.E_F).exp(); + damage_1 = 1. - ((self.parameters.E_0-kappa) / self.parameters.E_F).exp(); } else { damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); } diff --git a/src/jh2.rs b/src/jh2.rs index b90370b..975164f 100644 --- a/src/jh2.rs +++ b/src/jh2.rs @@ -26,7 +26,9 @@ pub struct JH2ConstParameters { pub EFMIN: f64, pub DMAX: f64, pub E_F: f64, + pub E_0: f64, pub LOCAL_SOUND_SPEED: f64, + pub HARDENING_MODULUS: f64, //pub REDUCE_T: f64, } #[derive(Debug)] @@ -58,7 +60,9 @@ impl ConstitutiveModel for JH23D { EFMIN: *parameters.get("EFMIN").unwrap(), DMAX: *parameters.get("DMAX").unwrap_or(&1.0), E_F: *parameters.get("E_F").unwrap_or(&0.0), + E_0: *parameters.get("E_0").unwrap_or(&0.0), LOCAL_SOUND_SPEED: *parameters.get("LOCAL_SOUND_SPEED").unwrap_or(&0.0), + HARDENING_MODULUS: *parameters.get("HARDENING_MODULUS").unwrap_or(&0.0), //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, }) @@ -118,7 +122,7 @@ impl ConstitutiveModel for JH23D { if self.parameters.E_F > 0.0 { let lambda_old = -self.parameters.E_F * (1. - damage_0).ln(); let lambda_new = lambda_old + del_lambda; - damage_1 = 1. - (-lambda_new / self.parameters.E_F).exp(); + damage_1 = 1. - ((self.parameters.E_0-lambda_new) / self.parameters.E_F).exp(); } else { damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); } From 40b4f04c1257039c149bb6d1cc7fa5455a309236 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 3 Apr 2024 16:02:24 +0200 Subject: [PATCH 20/24] Include hardening --- src/gradient_jh2.rs | 75 +++++++++++++++++++++++++-------------------- src/jh2.rs | 21 ++++++++----- 2 files changed, 54 insertions(+), 42 deletions(-) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index 40d9772..fc55a94 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -8,7 +8,7 @@ use crate::jh2::JH2ConstParameters; use std::collections::HashMap; #[derive(Debug)] -pub struct NonlocalParameters{ +pub struct NonlocalParameters { m_overnonlocal: f64, } @@ -19,7 +19,7 @@ pub struct GradientJH23D { } impl ConstitutiveModel for GradientJH23D { - fn new(parameters: &HashMap) -> Option{ + fn new(parameters: &HashMap) -> Option { Some(Self { parameters: JH2ConstParameters { RHO: *parameters.get("RHO").unwrap(), @@ -44,7 +44,7 @@ impl ConstitutiveModel for GradientJH23D { E_F: *parameters.get("E_F").unwrap_or(&0.0), E_0: *parameters.get("E_0").unwrap_or(&0.0), LOCAL_SOUND_SPEED: *parameters.get("LOCAL_SOUND_SPEED").unwrap_or(&0.0), - HARDENING_MODULUS: *parameters.get("HARDENING_MODULUS").unwrap_or(&0.0), + HARDENING: *parameters.get("HARDENING_MODULUS").unwrap_or(&0.0), //M_OVERNONLOCAL: *parameters.get("M_OVERNONLOCAL").unwrap_or(&0.0), //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, @@ -62,7 +62,7 @@ impl ConstitutiveModel for GradientJH23D { let d_eps = mandel_rate_from_velocity_gradient(&velocity_gradient); let (mut d_eps_vol, d_eps_dev) = mandel_decomposition(&d_eps); d_eps_vol *= -1.0; - + let sigma_0 = input.get_vector::<{ Q::MandelStress.size() }>(Q::MandelStress, ip); //let del_lambda_nonlocal = del_t * input.get_scalar(Q::EqNonlocalPlasticStrainRate, ip).max(0.0); let nonlocal_plastic_strain = input.get_scalar(Q::EqNonlocalPlasticStrain, ip); @@ -73,11 +73,10 @@ impl ConstitutiveModel for GradientJH23D { output.set_scalar(Q::HistoryMaximum, ip, history_maximum_1); let mut del_lambda = 0.0; - let (p_0, s_0) = mandel_decomposition(&sigma_0); let p_0 = p_0 - input.get_scalar(Q::BulkViscosity, ip); - + let s_tr = s_0 + 2. * self.parameters.SHEAR_MODULUS * d_eps_dev * del_t; let s_tr_eq = (1.5 * s_tr.norm_squared()).sqrt(); let d_eps_eq = ((2. / 3.) * d_eps.norm_squared()).sqrt(); @@ -86,26 +85,32 @@ impl ConstitutiveModel for GradientJH23D { let p_s = p_0 / self.parameters.PHEL; let t_s = self.parameters.T / self.parameters.PHEL; let mut rate_factor = 1.; - - - let e_p_f = (self.parameters.D1 * (p_s + t_s).powf(self.parameters.D2)) - .max(self.parameters.EFMIN); + + let e_p_f = + (self.parameters.D1 * (p_s + t_s).powf(self.parameters.D2)).max(self.parameters.EFMIN); let damage_0 = input.get_scalar(Q::Damage, ip); let mut damage_1 = 0.0; if self.parameters.E_F > 0.0 { let m = self.nonlocal_parameters.m_overnonlocal; - let kappa = (m * history_maximum_1 + (1.0 - m) * input.get_scalar(Q::EqPlasticStrain, ip)).max(0.0); - damage_1 = 1. - ((self.parameters.E_0-kappa) / self.parameters.E_F).exp(); + let kappa = (m * history_maximum_1 + + (1.0 - m) * input.get_scalar(Q::EqPlasticStrain, ip)) + .max(0.0); + damage_1 = 1. - ((self.parameters.E_0 - kappa).max(0.0) / self.parameters.E_F).exp(); } else { damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); } //let damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); output.set_scalar(Q::Damage, ip, damage_1); - + + let hardening = self.parameters.HARDENING; + let hardening_factor = hardening / ((1.0 - hardening) * self.parameters.E_0) + * input.get_scalar(Q::EqPlasticStrain, ip) + + 1.0; //let t_s_factor = (1.-damage_1).powf(self.parameters.REDUCE_T); + let yield_factor = 1.0 - hardening; let fracture_surface = - (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) + yield_factor * (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) .max(0.0); let residual_surface = (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); @@ -114,9 +119,9 @@ impl ConstitutiveModel for GradientJH23D { } let yield_surface = rate_factor * { if damage_0 == 0.0 { - fracture_surface + fracture_surface * hardening_factor } else { - fracture_surface * (1. - damage_1) + damage_1 * residual_surface + fracture_surface * hardening_factor * (1. - damage_1) + damage_1 * residual_surface } }; if s_tr_eq > yield_surface { @@ -125,8 +130,12 @@ impl ConstitutiveModel for GradientJH23D { } else { alpha = 1.0; } - - output.set_scalar(Q::EqPlasticStrain, ip, input.get_scalar(Q::EqPlasticStrain, ip) + del_lambda); + + output.set_scalar( + Q::EqPlasticStrain, + ip, + input.get_scalar(Q::EqPlasticStrain, ip) + del_lambda, + ); // /*********************************************************************** // * UPDATE DENSITY @@ -184,7 +193,7 @@ impl ConstitutiveModel for GradientJH23D { let c = (bulk_modulus / self.parameters.RHO).sqrt(); let q_1 = { if tr_deps < 0.0 { - l * density_1 * (1.5 * l * tr_deps.powi(2) + c * 0.06 * tr_deps.abs()) + l * density_1 * (1.5 * l * tr_deps.powi(2) + c * 0.06 * tr_deps.abs()) } else { 0.0 } @@ -195,7 +204,7 @@ impl ConstitutiveModel for GradientJH23D { // * Combine deviatoric and volumetric stresses // **********************************************************************/ let s_1 = s_tr * alpha; - let sigma_1 = s_1 - MANDEL_IDENTITY * (p_1+q_1); + let sigma_1 = s_1 - MANDEL_IDENTITY * (p_1 + q_1); output.set_vector(Q::MandelStress, ip, sigma_1); // *********************************************************************** @@ -214,39 +223,39 @@ impl ConstitutiveModel for GradientJH23D { if output.is_some(Q::Pressure) { output.set_scalar(Q::Pressure, ip, p_1); } - + let elastic_rate = - - (1. - alpha) / (2. * self.parameters.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; - + -(1. - alpha) / (2. * self.parameters.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; + let density_mid = 0.5 * (density_0 + density_1); if output.is_some(Q::InternalPlasticEnergy) && input.is_some(Q::InternalPlasticEnergy) { let s_mid = 0.5 * (s_0 + s_1); let deviatoric_rate = d_eps_dev - elastic_rate; let e_0 = input.get_scalar(Q::InternalPlasticEnergy, ip); - let e_1 = e_0 + del_t/density_mid * (s_mid.dot(&deviatoric_rate)); + let e_1 = e_0 + del_t / density_mid * (s_mid.dot(&deviatoric_rate)); output.set_scalar(Q::InternalPlasticEnergy, ip, e_1); } if output.is_some(Q::InternalElasticEnergy) && input.is_some(Q::InternalElasticEnergy) { let s_mid = 0.5 * (s_0 + s_1); - let p_mid = - 0.5 * (p_0 + p_1); + let p_mid = -0.5 * (p_0 + p_1); let deviatoric_rate = elastic_rate; let e_0 = input.get_scalar(Q::InternalElasticEnergy, ip); - let e_1 = e_0 + del_t/density_mid * (s_mid.dot(&deviatoric_rate) + 3. * d_eps_vol * p_mid); + let e_1 = + e_0 + del_t / density_mid * (s_mid.dot(&deviatoric_rate) + 3. * d_eps_vol * p_mid); output.set_scalar(Q::InternalElasticEnergy, ip, e_1); } if output.is_some(Q::InternalHeatingEnergy) && input.is_some(Q::InternalHeatingEnergy) { - let q_mid = - 0.5 * (q_1 + input.get_scalar(Q::BulkViscosity, ip)); + let q_mid = -0.5 * (q_1 + input.get_scalar(Q::BulkViscosity, ip)); let e_0 = input.get_scalar(Q::InternalHeatingEnergy, ip); - let e_1 = e_0 + del_t/density_mid * 3.* d_eps_vol * q_mid; + let e_1 = e_0 + del_t / density_mid * 3. * d_eps_vol * q_mid; output.set_scalar(Q::InternalHeatingEnergy, ip, e_1); } if output.is_some(Q::InternalEnergy) && input.is_some(Q::InternalEnergy) { let e_0 = input.get_scalar(Q::InternalEnergy, ip); let sigma_mid = 0.5 * (sigma_0 + sigma_1); - let e_1 = e_0 + del_t/density_mid * sigma_mid.dot(&d_eps); + let e_1 = e_0 + del_t / density_mid * sigma_mid.dot(&d_eps); output.set_scalar(Q::InternalEnergy, ip, e_1); } - } /// Returns the physical quantities that are required as input for the @@ -276,12 +285,10 @@ impl ConstitutiveModel for GradientJH23D { /// Returns the physical quantities that are needed as output, but are not /// necessarily needed in oredr to calculate the constitutive model. An example is - /// the consistent tangent which is not needed for the calculation of the stresses + /// the consistent tangent which is not needed for the calculation of the stresses /// and is therefore purely an output quantity. fn define_output(&self) -> HashMap { - HashMap::from([ - (Q::MandelStress, QDim::Vector(6)), - ]) + HashMap::from([(Q::MandelStress, QDim::Vector(6))]) } /// Returns the physical quantities that are optional output of the constitutive diff --git a/src/jh2.rs b/src/jh2.rs index 975164f..2961dd3 100644 --- a/src/jh2.rs +++ b/src/jh2.rs @@ -28,7 +28,7 @@ pub struct JH2ConstParameters { pub E_F: f64, pub E_0: f64, pub LOCAL_SOUND_SPEED: f64, - pub HARDENING_MODULUS: f64, + pub HARDENING: f64, //pub REDUCE_T: f64, } #[derive(Debug)] @@ -62,7 +62,7 @@ impl ConstitutiveModel for JH23D { E_F: *parameters.get("E_F").unwrap_or(&0.0), E_0: *parameters.get("E_0").unwrap_or(&0.0), LOCAL_SOUND_SPEED: *parameters.get("LOCAL_SOUND_SPEED").unwrap_or(&0.0), - HARDENING_MODULUS: *parameters.get("HARDENING_MODULUS").unwrap_or(&0.0), + HARDENING: *parameters.get("HARDENING").unwrap_or(&0.0), //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, }) @@ -83,6 +83,10 @@ impl ConstitutiveModel for JH23D { let damage_0 = input.get_scalar(Q::Damage, ip); let mut damage_1 = damage_0; + + let lambda_old = -self.parameters.E_F * (1. - damage_0).ln(); + let hardening = self.parameters.HARDENING; + let hardening_factor = hardening/((1.0-hardening) * self.parameters.E_0) * lambda_old + 1.0; let (p_0, s_0) = mandel_decomposition(&sigma_0); let p_0 = p_0 - input.get_scalar(Q::BulkViscosity, ip); @@ -90,15 +94,16 @@ impl ConstitutiveModel for JH23D { let s_tr = s_0 + 2. * self.parameters.SHEAR_MODULUS * d_eps_dev * del_t; let s_tr_eq = (1.5 * s_tr.norm_squared()).sqrt(); let d_eps_eq = ((2. / 3.) * d_eps.norm_squared()).sqrt(); - let mut alpha; + let mut alpha: f64; let p_s = p_0 / self.parameters.PHEL; let t_s = self.parameters.T / self.parameters.PHEL; let mut rate_factor = 1.; //let t_s_factor = (1.-damage_0).powf(self.parameters.REDUCE_T); + let yield_factor = 1.0 - hardening; let fracture_surface = - (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) + yield_factor * (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) .max(0.0); let residual_surface = (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); @@ -107,9 +112,9 @@ impl ConstitutiveModel for JH23D { } let yield_surface = rate_factor * { if damage_0 == 0.0 { - fracture_surface + fracture_surface * hardening_factor } else { - fracture_surface * (1. - damage_0) + damage_0 * residual_surface + fracture_surface * hardening_factor * (1. - damage_0) + damage_0 * residual_surface } }; if s_tr_eq > yield_surface { @@ -120,9 +125,9 @@ impl ConstitutiveModel for JH23D { alpha = yield_surface / s_tr_eq; if self.parameters.E_F > 0.0 { - let lambda_old = -self.parameters.E_F * (1. - damage_0).ln(); + //let lambda_old = -self.parameters.E_F * (1. - damage_0).ln(); let lambda_new = lambda_old + del_lambda; - damage_1 = 1. - ((self.parameters.E_0-lambda_new) / self.parameters.E_F).exp(); + damage_1 = 1. - ((self.parameters.E_0-lambda_new).max(0.0) / self.parameters.E_F).exp(); } else { damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); } From 755c3fa1bb80f1ec2c8f337215c4878de94da808 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 3 Apr 2024 20:01:13 +0200 Subject: [PATCH 21/24] Fix nan error --- src/gradient_jh2.rs | 13 +++++++++---- src/jh2.rs | 20 +++++++++++++++----- 2 files changed, 24 insertions(+), 9 deletions(-) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index fc55a94..dc13d3b 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -44,7 +44,7 @@ impl ConstitutiveModel for GradientJH23D { E_F: *parameters.get("E_F").unwrap_or(&0.0), E_0: *parameters.get("E_0").unwrap_or(&0.0), LOCAL_SOUND_SPEED: *parameters.get("LOCAL_SOUND_SPEED").unwrap_or(&0.0), - HARDENING: *parameters.get("HARDENING_MODULUS").unwrap_or(&0.0), + HARDENING: *parameters.get("HARDENING").unwrap_or(&0.0), //M_OVERNONLOCAL: *parameters.get("M_OVERNONLOCAL").unwrap_or(&0.0), //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, @@ -104,9 +104,14 @@ impl ConstitutiveModel for GradientJH23D { output.set_scalar(Q::Damage, ip, damage_1); let hardening = self.parameters.HARDENING; - let hardening_factor = hardening / ((1.0 - hardening) * self.parameters.E_0) - * input.get_scalar(Q::EqPlasticStrain, ip) - + 1.0; + let hardening_slope = hardening / ((1.0 - hardening) * self.parameters.E_0); + let hardening_factor = { + if hardening_slope.is_nan() { + 1.0 + } else { + hardening_slope * input.get_scalar(Q::EqPlasticStrain, ip) + 1.0 + } + }; //let t_s_factor = (1.-damage_1).powf(self.parameters.REDUCE_T); let yield_factor = 1.0 - hardening; let fracture_surface = diff --git a/src/jh2.rs b/src/jh2.rs index 2961dd3..ce0a467 100644 --- a/src/jh2.rs +++ b/src/jh2.rs @@ -3,6 +3,7 @@ use crate::stress_strain::{ mandel_decomposition, mandel_rate_from_velocity_gradient, MANDEL_IDENTITY, }; +use core::panic; use std::collections::HashMap; #[derive(Debug)] pub struct JH2ConstParameters { @@ -83,10 +84,18 @@ impl ConstitutiveModel for JH23D { let damage_0 = input.get_scalar(Q::Damage, ip); let mut damage_1 = damage_0; - + let lambda_old = -self.parameters.E_F * (1. - damage_0).ln(); let hardening = self.parameters.HARDENING; - let hardening_factor = hardening/((1.0-hardening) * self.parameters.E_0) * lambda_old + 1.0; + + let hardening_slope = hardening / ((1.0 - hardening) * self.parameters.E_0); + let hardening_factor = { + if hardening_slope.is_nan() { + 1.0 + } else { + hardening_slope * lambda_old + 1.0 + } + }; let (p_0, s_0) = mandel_decomposition(&sigma_0); let p_0 = p_0 - input.get_scalar(Q::BulkViscosity, ip); @@ -102,8 +111,8 @@ impl ConstitutiveModel for JH23D { //let t_s_factor = (1.-damage_0).powf(self.parameters.REDUCE_T); let yield_factor = 1.0 - hardening; - let fracture_surface = - yield_factor * (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) + let fracture_surface = yield_factor + * (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) .max(0.0); let residual_surface = (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); @@ -127,7 +136,8 @@ impl ConstitutiveModel for JH23D { if self.parameters.E_F > 0.0 { //let lambda_old = -self.parameters.E_F * (1. - damage_0).ln(); let lambda_new = lambda_old + del_lambda; - damage_1 = 1. - ((self.parameters.E_0-lambda_new).max(0.0) / self.parameters.E_F).exp(); + damage_1 = + 1. - ((self.parameters.E_0 - lambda_new).max(0.0) / self.parameters.E_F).exp(); } else { damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); } From dc43fd37e8d0858acce32560abedfc7343a141d9 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 3 Apr 2024 20:17:05 +0200 Subject: [PATCH 22/24] logic error in damage --- src/gradient_jh2.rs | 2 +- src/jh2.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index dc13d3b..9f3bae6 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -96,7 +96,7 @@ impl ConstitutiveModel for GradientJH23D { let kappa = (m * history_maximum_1 + (1.0 - m) * input.get_scalar(Q::EqPlasticStrain, ip)) .max(0.0); - damage_1 = 1. - ((self.parameters.E_0 - kappa).max(0.0) / self.parameters.E_F).exp(); + damage_1 = (1. - ((self.parameters.E_0 - kappa) / self.parameters.E_F).exp()).max(0.0); } else { damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); } diff --git a/src/jh2.rs b/src/jh2.rs index ce0a467..6563f27 100644 --- a/src/jh2.rs +++ b/src/jh2.rs @@ -137,7 +137,7 @@ impl ConstitutiveModel for JH23D { //let lambda_old = -self.parameters.E_F * (1. - damage_0).ln(); let lambda_new = lambda_old + del_lambda; damage_1 = - 1. - ((self.parameters.E_0 - lambda_new).max(0.0) / self.parameters.E_F).exp(); + (1. - ((self.parameters.E_0 - lambda_new) / self.parameters.E_F).exp()).max(0.0); } else { damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); } From 0d84a8f8c2bc3f0ed2b6f0ddc431f002b2c0a8be Mon Sep 17 00:00:00 2001 From: srosenbu Date: Thu, 4 Apr 2024 13:01:12 +0200 Subject: [PATCH 23/24] Fix lambda calculation --- src/jh2.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jh2.rs b/src/jh2.rs index 6563f27..6825ce8 100644 --- a/src/jh2.rs +++ b/src/jh2.rs @@ -85,7 +85,7 @@ impl ConstitutiveModel for JH23D { let damage_0 = input.get_scalar(Q::Damage, ip); let mut damage_1 = damage_0; - let lambda_old = -self.parameters.E_F * (1. - damage_0).ln(); + let lambda_old = -self.parameters.E_F * (1. - damage_0).ln() + self.parameters.E_0; let hardening = self.parameters.HARDENING; let hardening_slope = hardening / ((1.0 - hardening) * self.parameters.E_0); From da66e2bb26f98aea86a422464b5c7210fe3e34fb Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 14 Aug 2024 09:57:03 +0200 Subject: [PATCH 24/24] Add Cargo lock file to version control --- Cargo.lock | 496 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 496 insertions(+) create mode 100644 Cargo.lock diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..7065a30 --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,496 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + +[[package]] +name = "autocfg" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" + +[[package]] +name = "bitflags" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" + +[[package]] +name = "bytemuck" +version = "1.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "17febce684fd15d89027105661fec94afb475cb995fbc59d2865198446ba2eea" + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "comfe" +version = "0.1.0" +dependencies = [ + "nalgebra", + "numpy", + "pyo3", + "simba", + "strum", + "strum_macros", +] + +[[package]] +name = "heck" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" + +[[package]] +name = "indoc" +version = "1.0.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bfa799dd5ed20a7e349f3b4639aa80d74549c81716d9ec4f994c9b5815598306" + +[[package]] +name = "libc" +version = "0.2.147" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b4668fb0ea861c1df094127ac5f1da3409a82116a4ba74fca2e58ef927159bb3" + +[[package]] +name = "lock_api" +version = "0.4.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c1cc9717a20b1bb222f333e6a92fd32f7d8a18ddc5a3191a11af45dcbf4dcd16" +dependencies = [ + "autocfg", + "scopeguard", +] + +[[package]] +name = "matrixmultiply" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "090126dc04f95dc0d1c1c91f61bdd474b3930ca064c1edc8a849da2c6cbe1e77" +dependencies = [ + "autocfg", + "rawpointer", +] + +[[package]] +name = "memoffset" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a634b1c61a95585bd15607c6ab0c4e5b226e695ff2800ba0cdccddf208c406c" +dependencies = [ + "autocfg", +] + +[[package]] +name = "nalgebra" +version = "0.32.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "307ed9b18cc2423f29e83f84fd23a8e73628727990181f18641a8b5dc2ab1caa" +dependencies = [ + "approx", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "nalgebra-macros" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91761aed67d03ad966ef783ae962ef9bbaca728d2dd7ceb7939ec110fffad998" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "ndarray" +version = "0.15.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adb12d4e967ec485a5f71c6311fe28158e9d6f4bc4a447b474184d0f91a8fa32" +dependencies = [ + "matrixmultiply", + "num-complex", + "num-integer", + "num-traits", + "rawpointer", +] + +[[package]] +name = "num-complex" +version = "0.4.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ba157ca0885411de85d6ca030ba7e2a83a28636056c7c699b07c8b6f7383214" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" +dependencies = [ + "autocfg", + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f30b0abd723be7e2ffca1272140fac1a2f084c77ec3e123c192b66af1ee9e6c2" +dependencies = [ + "autocfg", +] + +[[package]] +name = "numpy" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "437213adf41bbccf4aeae535fbfcdad0f6fed241e1ae182ebe97fa1f3ce19389" +dependencies = [ + "libc", + "nalgebra", + "ndarray", + "num-complex", + "num-integer", + "num-traits", + "pyo3", + "rustc-hash", +] + +[[package]] +name = "once_cell" +version = "1.18.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dd8b5dd2ae5ed71462c540258bedcb51965123ad7e7ccf4b9a8cafaa4a63576d" + +[[package]] +name = "parking_lot" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3742b2c103b9f06bc9fff0a37ff4912935851bee6d36f3c02bcc755bcfec228f" +dependencies = [ + "lock_api", + "parking_lot_core", +] + +[[package]] +name = "parking_lot_core" +version = "0.9.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93f00c865fe7cabf650081affecd3871070f26767e7b2070a3ffae14c654b447" +dependencies = [ + "cfg-if", + "libc", + "redox_syscall", + "smallvec", + "windows-targets", +] + +[[package]] +name = "paste" +version = "1.0.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "de3145af08024dea9fa9914f381a17b8fc6034dfb00f3a84013f7ff43f29ed4c" + +[[package]] +name = "proc-macro2" +version = "1.0.66" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "18fb31db3f9bddb2ea821cde30a9f70117e3f119938b5ee630b7403aa6e2ead9" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pyo3" +version = "0.19.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e681a6cfdc4adcc93b4d3cf993749a4552018ee0a9b65fc0ccfad74352c72a38" +dependencies = [ + "cfg-if", + "indoc", + "libc", + "memoffset", + "parking_lot", + "pyo3-build-config", + "pyo3-ffi", + "pyo3-macros", + "unindent", +] + +[[package]] +name = "pyo3-build-config" +version = "0.19.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "076c73d0bc438f7a4ef6fdd0c3bb4732149136abd952b110ac93e4edb13a6ba5" +dependencies = [ + "once_cell", + "target-lexicon", +] + +[[package]] +name = "pyo3-ffi" +version = "0.19.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e53cee42e77ebe256066ba8aa77eff722b3bb91f3419177cf4cd0f304d3284d9" +dependencies = [ + "libc", + "pyo3-build-config", +] + +[[package]] +name = "pyo3-macros" +version = "0.19.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dfeb4c99597e136528c6dd7d5e3de5434d1ceaf487436a3f03b2d56b6fc9efd1" +dependencies = [ + "proc-macro2", + "pyo3-macros-backend", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "pyo3-macros-backend" +version = "0.19.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "947dc12175c254889edc0c02e399476c2f652b4b9ebd123aa655c224de259536" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "quote" +version = "1.0.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5267fca4496028628a95160fc423a33e8b2e6af8a5302579e322e4b520293cae" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + +[[package]] +name = "redox_syscall" +version = "0.3.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "567664f262709473930a4bf9e51bf2ebf3348f2e748ccc50dea20646858f8f29" +dependencies = [ + "bitflags", +] + +[[package]] +name = "rustc-hash" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" + +[[package]] +name = "rustversion" +version = "1.0.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ffc183a10b4478d04cbbbfc96d0873219d962dd5accaff2ffbd4ceb7df837f4" + +[[package]] +name = "safe_arch" +version = "0.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f398075ce1e6a179b46f51bd88d0598b92b00d3551f1a2d4ac49e771b56ac354" +dependencies = [ + "bytemuck", +] + +[[package]] +name = "scopeguard" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" + +[[package]] +name = "simba" +version = "0.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "061507c94fc6ab4ba1c9a0305018408e312e17c041eb63bef8aa726fa33aceae" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "paste", + "wide", +] + +[[package]] +name = "smallvec" +version = "1.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62bb4feee49fdd9f707ef802e22365a35de4b7b299de4763d44bfea899442ff9" + +[[package]] +name = "strum" +version = "0.25.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "290d54ea6f91c969195bdbcd7442c8c2a2ba87da8bf60a7ee86a235d4bc1e125" + +[[package]] +name = "strum_macros" +version = "0.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ad8d03b598d3d0fff69bf533ee3ef19b8eeb342729596df84bcc7e1f96ec4059" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "rustversion", + "syn 2.0.29", +] + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "2.0.29" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c324c494eba9d92503e6f1ef2e6df781e78f6a7705a0202d9801b198807d518a" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "target-lexicon" +version = "0.12.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9d0e916b1148c8e263850e1ebcbd046f333e0683c724876bb0da63ea4373dc8a" + +[[package]] +name = "typenum" +version = "1.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "497961ef93d974e23eb6f433eb5fe1b7930b659f06d12dec6fc44a8f554c0bba" + +[[package]] +name = "unicode-ident" +version = "1.0.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "301abaae475aa91687eb82514b328ab47a211a533026cb25fc3e519b86adfc3c" + +[[package]] +name = "unindent" +version = "0.1.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e1766d682d402817b5ac4490b3c3002d91dfa0d22812f341609f97b08757359c" + +[[package]] +name = "wide" +version = "0.7.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa469ffa65ef7e0ba0f164183697b89b854253fd31aeb92358b7b6155177d62f" +dependencies = [ + "bytemuck", + "safe_arch", +] + +[[package]] +name = "windows-targets" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9a2fa6e2155d7247be68c096456083145c183cbbbc2764150dda45a87197940c" +dependencies = [ + "windows_aarch64_gnullvm", + "windows_aarch64_msvc", + "windows_i686_gnu", + "windows_i686_msvc", + "windows_x86_64_gnu", + "windows_x86_64_gnullvm", + "windows_x86_64_msvc", +] + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b38e32f0abccf9987a4e3079dfb67dcd799fb61361e53e2882c3cbaf0d905d8" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dc35310971f3b2dbbf3f0690a219f40e2d9afcf64f9ab7cc1be722937c26b4bc" + +[[package]] +name = "windows_i686_gnu" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a75915e7def60c94dcef72200b9a8e58e5091744960da64ec734a6c6e9b3743e" + +[[package]] +name = "windows_i686_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f55c233f70c4b27f66c523580f78f1004e8b5a8b659e05a4eb49d4166cca406" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "53d40abd2583d23e4718fddf1ebec84dbff8381c07cae67ff7768bbf19c6718e" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b7b52767868a23d5bab768e390dc5f5c55825b6d30b86c844ff2dc7414044cc" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed94fce61571a4006852b7389a063ab983c02eb1bb37b47f8272ce92d06d9538"