Detect invalid number

TLCFEM · Oct 17, 2023 · 53e94be · 53e94be
1 parent 0b14d1e
commit 53e94be
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 12 deletions.
diff --git a/Material/Material3D/Concrete/CDPM2.cpp b/Material/Material3D/Concrete/CDPM2.cpp
@@ -25,7 +25,7 @@ const double CDPM2::sqrt_six = std::sqrt(6.);
 const double CDPM2::sqrt_three_two = std::sqrt(1.5);
 const mat CDPM2::unit_dev_tensor = tensor::unit_deviatoric_tensor4();
 
-void CDPM2::compute_plasticity(const double lode, const double s, const double p, const double kp, podarray<double>& data) const {
+void CDPM2::compute_plasticity(const double lode, const double s, const double p, const double kp, vec& data) const {
     auto& f = data(0);
     auto& pfps = data(1);
     auto& pfpp = data(2);
@@ -131,9 +131,11 @@ void CDPM2::compute_plasticity(const double lode, const double s, const double p
         dxhdp = -xh / fh / fc;
         xh += dh;
     }
+
+    suanpan_assert([&] { if(!data.is_finite()) throw std::runtime_error("unexpected non-finite internal data"); });
 }
 
-int CDPM2::compute_damage(const double gamma, const double s, const double p, const double kp, const double ac, podarray<double>& data) {
+int CDPM2::compute_damage(const double gamma, const double s, const double p, const double kp, const double ac, vec& data) {
     const auto& gs = data(4);
     const auto& gp = data(5);
     const auto& gg = data(6);
@@ -255,7 +257,7 @@ int CDPM2::compute_damage(const double gamma, const double s, const double p, co
     const auto pkdc2pp = (pkdcpp - incre_kdc2 * pxspp) / xs;
     const auto pkdc2pac = pkdcpac / xs;
 
-    podarray<double> datad(3);
+    vec datad(3);
 
     if(SUANPAN_SUCCESS != compute_damage_factor(kdt, kdt1, kdt2, eft, omegat, datad)) return SUANPAN_FAIL;
     const auto& potpkdt = datad(0);
@@ -292,7 +294,7 @@ int CDPM2::compute_damage(const double gamma, const double s, const double p, co
     return SUANPAN_SUCCESS;
 }
 
-int CDPM2::compute_damage_factor(const double kd, const double kd1, const double kd2, const double ef, double& omega, podarray<double>& data) const {
+int CDPM2::compute_damage_factor(const double kd, const double kd1, const double kd2, const double ef, double& omega, vec& data) const {
     auto& popkd = data(0);
     auto& popkd1 = data(1);
     auto& popkd2 = data(2);
@@ -394,7 +396,7 @@ int CDPM2::update_trial_status(const vec& t_strain) {
 
     vec residual(4), incre;
 
-    podarray<double> data(16);
+    vec data(16);
     const auto& f = data(0);
     const auto& pfps = data(1);
     const auto& pfpp = data(2);
@@ -429,6 +431,8 @@ int CDPM2::update_trial_status(const vec& t_strain) {
                     return f;
                 };
 
+                approx_update(0.); // clear data
+
                 gamma = ini_f / elastic_modulus;
                 auto x1 = 0., f1 = ini_f;
                 // find a proper bracket
@@ -439,13 +443,13 @@ int CDPM2::update_trial_status(const vec& t_strain) {
                     gamma *= 2.;
                 }
 
-                auto t_counter = 0u;
+                counter = 0u;
                 auto x2 = gamma, f2 = f;
                 while(true) {
                     // Ridders' method
                     if(fabs(x2 - x1) < tolerance) break;
 
-                    t_counter += 2u;
+                    counter += 2u;
 
                     const auto x3 = .5 * (x1 + x2);
                     const auto f3 = approx_update(x3);
@@ -472,9 +476,9 @@ int CDPM2::update_trial_status(const vec& t_strain) {
                     f2 = f4;
                 }
 
-                counter = 1u; // avoid initial elastic check
+                suanpan_debug("Ridders' method initial guess {:.5E} with {} iterations.\n", gamma, counter);
 
-                suanpan_debug("Ridders' method initial guess {:.10E} with {} iterations.\n", gamma, t_counter);
+                counter = 1u; // avoid initial elastic check
 
                 continue;
             }

diff --git a/Material/Material3D/Concrete/CDPM2.h b/Material/Material3D/Concrete/CDPM2.h
@@ -86,9 +86,9 @@ class CDPM2 final : protected DataCDPM2, public Material3D {
 
     const DamageType damage_type = DamageType::ANISOTROPIC;
 
-    void compute_plasticity(double, double, double, double, podarray<double>&) const;
-    int compute_damage(double, double, double, double, double, podarray<double>&);
-    int compute_damage_factor(double, double, double, double, double&, podarray<double>&) const;
+    void compute_plasticity(double, double, double, double, vec&) const;
+    int compute_damage(double, double, double, double, double, vec&);
+    int compute_damage_factor(double, double, double, double, double&, vec&) const;
 
 public:
     CDPM2(unsigned,   // tag