diff --git a/README.md b/README.md index 0cb6b345..194921f2 100644 --- a/README.md +++ b/README.md @@ -42,6 +42,7 @@ conveniently compiled to this [PDF](docs/readme.pdf). - [Upgrading instructions](#upgrading-instructions) - [Upgrading to v0.6.*](#upgrading-to-v06) - [Change-log](#change-log) + - [v0.7.0](#v070) - [v0.6.4](#v064) - [v0.6.3](#v063) - [v0.6.2](#v062) @@ -466,6 +467,11 @@ This requires the following changes: # Change-log +## v0.7.0 + +* Exposing "checkYieldBoundLeft" and "checkYieldBoundRight" from QPot v0.2.0. +* Proper float testing using Catch2. + ## v0.6.4 * Increasing robustness using the new xtensor's `has_fixed_rank` and `get_rank`. diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index cea96bc8..a8962d70 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -21,6 +21,8 @@ foreach(testsourcefile ${APP_SOURCES}) add_executable(${testname} ${testsourcefile}) find_package(GooseFEM REQUIRED) + find_package(xtensor REQUIRED) + find_package(xsimd REQUIRED) target_link_libraries(${testname} PUBLIC GMatElastoPlasticQPot diff --git a/examples/friction-optimized.cpp b/examples/friction-optimized.cpp index bcdaa2e1..dc46f33b 100644 --- a/examples/friction-optimized.cpp +++ b/examples/friction-optimized.cpp @@ -314,6 +314,8 @@ class System { xt::noalias(m_u) = m_u + m_dt * m_v + 0.5 * std::pow(m_dt, 2.) * m_a; + // compute strain/strain, and corresponding force + computeStrainStressForcesWeakLayer(); #ifdef CHECK_MYSHORTCUT diff --git a/examples/friction.cpp b/examples/friction.cpp index f3c99965..6d9cb1e3 100644 --- a/examples/friction.cpp +++ b/examples/friction.cpp @@ -140,6 +140,7 @@ class System { m_Eps = m_quad.AllocateQtensor<2>(0.0); m_Sig = m_quad.AllocateQtensor<2>(0.0); + // -------- // material // -------- diff --git a/include/GMatElastoPlasticQPot/Cartesian2d.h b/include/GMatElastoPlasticQPot/Cartesian2d.h index 73eda842..a21f6050 100644 --- a/include/GMatElastoPlasticQPot/Cartesian2d.h +++ b/include/GMatElastoPlasticQPot/Cartesian2d.h @@ -149,6 +149,10 @@ class Cusp double epsp() const; // "plastic strain" (mean of currentYieldLeft and currentYieldRight) double energy() const; // potential energy + // Check that 'the particle' is at least "n" wells from the far-left/right + bool checkYieldBoundLeft(size_t n = 0) const; + bool checkYieldBoundRight(size_t n = 0) const; + private: double m_K; // bulk modulus double m_G; // shear modulus @@ -200,6 +204,10 @@ class Smooth double epsp() const; // "plastic strain" (mean of currentYieldLeft and currentYieldRight) double energy() const; // potential energy + // Check that 'the particle' is at least "n" wells from the far-left/right + bool checkYieldBoundLeft(size_t n = 0) const; + bool checkYieldBoundRight(size_t n = 0) const; + private: double m_K; // bulk modulus double m_G; // shear modulus @@ -316,6 +324,8 @@ class Array void currentIndex(xt::xtensor& arg) const; void currentYieldLeft(xt::xtensor& arg) const; void currentYieldRight(xt::xtensor& arg) const; + bool checkYieldBoundLeft(size_t n = 0) const; + bool checkYieldBoundRight(size_t n = 0) const; void epsp(xt::xtensor& arg) const; void energy(xt::xtensor& arg) const; diff --git a/include/GMatElastoPlasticQPot/Cartesian2d_Array.hpp b/include/GMatElastoPlasticQPot/Cartesian2d_Array.hpp index 64112a74..ee631c48 100644 --- a/include/GMatElastoPlasticQPot/Cartesian2d_Array.hpp +++ b/include/GMatElastoPlasticQPot/Cartesian2d_Array.hpp @@ -502,6 +502,58 @@ inline void Array::currentIndex(xt::xtensor& A) const } } +template +inline bool Array::checkYieldBoundLeft(size_t n) const +{ + GMATELASTOPLASTICQPOT_ASSERT(m_allSet); + + #pragma omp parallel for + for (size_t i = 0; i < m_size; ++i) { + switch (m_type.data()[i]) { + case Type::Elastic: + break; + case Type::Cusp: + if (!m_Cusp[m_index.data()[i]].checkYieldBoundLeft(n)) { + return false; + } + break; + case Type::Smooth: + if (!m_Smooth[m_index.data()[i]].checkYieldBoundLeft(n)) { + return false; + } + break; + } + } + + return true; +} + +template +inline bool Array::checkYieldBoundRight(size_t n) const +{ + GMATELASTOPLASTICQPOT_ASSERT(m_allSet); + + #pragma omp parallel for + for (size_t i = 0; i < m_size; ++i) { + switch (m_type.data()[i]) { + case Type::Elastic: + break; + case Type::Cusp: + if (!m_Cusp[m_index.data()[i]].checkYieldBoundRight(n)) { + return false; + } + break; + case Type::Smooth: + if (!m_Smooth[m_index.data()[i]].checkYieldBoundRight(n)) { + return false; + } + break; + } + } + + return true; +} + template inline void Array::currentYieldLeft(xt::xtensor& A) const { diff --git a/include/GMatElastoPlasticQPot/Cartesian2d_Cusp.hpp b/include/GMatElastoPlasticQPot/Cartesian2d_Cusp.hpp index f806b3d4..c43106bb 100644 --- a/include/GMatElastoPlasticQPot/Cartesian2d_Cusp.hpp +++ b/include/GMatElastoPlasticQPot/Cartesian2d_Cusp.hpp @@ -147,6 +147,16 @@ inline double Cusp::energy() const return U + V; } +inline bool Cusp::checkYieldBoundLeft(size_t n) const +{ + return m_yield.checkYieldBoundLeft(n); +} + +inline bool Cusp::checkYieldBoundRight(size_t n) const +{ + return m_yield.checkYieldBoundRight(n); +} + } // namespace Cartesian2d } // namespace GMatElastoPlasticQPot diff --git a/include/GMatElastoPlasticQPot/Cartesian2d_Smooth.hpp b/include/GMatElastoPlasticQPot/Cartesian2d_Smooth.hpp index 04448a3f..318c9f96 100644 --- a/include/GMatElastoPlasticQPot/Cartesian2d_Smooth.hpp +++ b/include/GMatElastoPlasticQPot/Cartesian2d_Smooth.hpp @@ -149,6 +149,16 @@ inline double Smooth::energy() const return U + V; } +inline bool Smooth::checkYieldBoundLeft(size_t n) const +{ + return m_yield.checkYieldBoundLeft(n); +} + +inline bool Smooth::checkYieldBoundRight(size_t n) const +{ + return m_yield.checkYieldBoundRight(n); +} + } // namespace Cartesian2d } // namespace GMatElastoPlasticQPot diff --git a/include/GMatElastoPlasticQPot/config.h b/include/GMatElastoPlasticQPot/config.h index f5abf4dc..5f014cd8 100644 --- a/include/GMatElastoPlasticQPot/config.h +++ b/include/GMatElastoPlasticQPot/config.h @@ -49,8 +49,8 @@ #endif #define GMATELASTOPLASTICQPOT_VERSION_MAJOR 0 -#define GMATELASTOPLASTICQPOT_VERSION_MINOR 6 -#define GMATELASTOPLASTICQPOT_VERSION_PATCH 4 +#define GMATELASTOPLASTICQPOT_VERSION_MINOR 7 +#define GMATELASTOPLASTICQPOT_VERSION_PATCH 0 #define GMATELASTOPLASTICQPOT_VERSION_AT_LEAST(x,y,z) \ (GMATELASTOPLASTICQPOT_VERSION_MAJOR > x || (GMATELASTOPLASTICQPOT_VERSION_MAJOR >= x && \ diff --git a/python/main.cpp b/python/main.cpp index 508f13bb..3da76ed2 100644 --- a/python/main.cpp +++ b/python/main.cpp @@ -24,59 +24,35 @@ auto add_common_members_array(T& self) self.def(py::init>(), "Matrix of material points.", py::arg("shape")) - .def("shape", &SM::Array::shape, "Return matrix shape.") + .def("shape", &SM::Array::shape, "Matrix shape.") - .def("K", &SM::Array::K, "Return matrix with bulk moduli.") + .def("K", &SM::Array::K, "Matrix with bulk moduli.") - .def("G", &SM::Array::G, "Return matrix with shear moduli.") + .def("G", &SM::Array::G, "Matrix with shear moduli.") - .def("I2", &SM::Array::I2, "Return matrix with second order unit tensors.") + .def("I2", &SM::Array::I2, "Matrix with 2nd-order unit tensors.") - .def( - "II", - &SM::Array::II, - "Return matrix with fourth order tensors with the result of the dyadic product II.") + .def("II", &SM::Array::II, "Matrix with 4th-order tensors = dyadic(I2, I2).") - .def("I4", &SM::Array::I4, "Return matrix with fourth order unit tensors.") + .def("I4", &SM::Array::I4, "Matrix with 4th-order unit tensors.") - .def( - "I4rt", - &SM::Array::I4rt, - "Return matrix with fourth right-transposed order unit tensors.") + .def("I4rt", &SM::Array::I4rt, "Matrix with 4th-order right-transposed unit tensors.") - .def( - "I4s", - &SM::Array::I4s, - "Return matrix with fourth order symmetric projection tensors.") + .def("I4s", &SM::Array::I4s, "Matrix with 4th-order symmetric projection tensors.") - .def( - "I4d", - &SM::Array::I4d, - "Return matrix with fourth order deviatoric projection tensors.") + .def("I4d", &SM::Array::I4d, "Matrix with 4th-order deviatoric projection tensors.") - .def("type", &SM::Array::type, "Return matrix with material types.") + .def("type", &SM::Array::type, "Matrix with material types.") - .def( - "isElastic", - &SM::Array::isElastic, - "Return matrix with boolean: Elastic (1) or not (0).") + .def("isElastic", &SM::Array::isElastic, "Boolean-matrix: true for Elastic.") - .def( - "isPlastic", - &SM::Array::isPlastic, - "Return matrix with boolean: Elastic (0) or plastic (Cusp/Smooth) (1).") + .def("isPlastic", &SM::Array::isPlastic, "Boolean-matrix: true for Cusp/Smooth.") - .def("isCusp", &SM::Array::isCusp, "Return matrix with boolean: Cusp (1) or not (0).") + .def("isCusp", &SM::Array::isCusp, "Boolean-matrix: true for Cusp.") - .def( - "isSmooth", - &SM::Array::isSmooth, - "Return matrix with boolean: Smooth (1) or not (0).") + .def("isSmooth", &SM::Array::isSmooth, "Boolean-matrix: true for Smooth.") - .def( - "check", - &SM::Array::check, - "Check that all matrix entries are set. Throws if any unset point is found.") + .def("check", &SM::Array::check, "Throws if any unset point is found.") .def( "setElastic", @@ -195,6 +171,18 @@ auto add_common_members_array(T& self) &SM::Array::CurrentYieldRight, "Returns matrix of yield strains to the left, given the current strain.") + .def( + "checkYieldBoundLeft", + &SM::Array::checkYieldBoundLeft, + "Check that 'the particle' is at least 'n' wells from the far-left.", + py::arg("n") = 0) + + .def( + "checkYieldBoundRight", + &SM::Array::checkYieldBoundRight, + "Check that 'the particle' is at least 'n' wells from the far-right.", + py::arg("n") = 0) + .def( "Epsp", &SM::Array::Epsp, @@ -385,6 +373,18 @@ PYBIND11_MODULE(GMatElastoPlasticQPot, m) &SM::Cusp::currentYieldRight, "Returns the yield strain to the right, for last known strain.") + .def( + "checkYieldBoundLeft", + &SM::Cusp::checkYieldBoundLeft, + "Check that 'the particle' is at least 'n' wells from the far-left.", + py::arg("n") = 0) + + .def( + "checkYieldBoundRight", + &SM::Cusp::checkYieldBoundRight, + "Check that 'the particle' is at least 'n' wells from the far-right.", + py::arg("n") = 0) + .def( "epsp", &SM::Cusp::epsp, "Returns the equivalent plastic strain for last known strain.") @@ -436,6 +436,18 @@ PYBIND11_MODULE(GMatElastoPlasticQPot, m) &SM::Smooth::currentYieldRight, "Returns the yield strain to the right, for last known strain.") + .def( + "checkYieldBoundLeft", + &SM::Smooth::checkYieldBoundLeft, + "Check that 'the particle' is at least 'n' wells from the far-left.", + py::arg("n") = 0) + + .def( + "checkYieldBoundRight", + &SM::Smooth::checkYieldBoundRight, + "Check that 'the particle' is at least 'n' wells from the far-right.", + py::arg("n") = 0) + .def( "epsp", &SM::Smooth::epsp, diff --git a/test/Cartesian2d.cpp b/test/Cartesian2d.cpp index a0d66b24..17476e00 100644 --- a/test/Cartesian2d.cpp +++ b/test/Cartesian2d.cpp @@ -1,11 +1,10 @@ #include #include +#include #define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); -#include - namespace GM = GMatElastoPlasticQPot::Cartesian2d; template @@ -89,7 +88,7 @@ SECTION("Hydrostatic - Tensor2") GM::Tensor2 A = xt::random::randn({2, 2}); A(0, 0) = 1.0; A(1, 1) = 1.0; - ISCLOSE(GM::Hydrostatic(A)(), 1.0); + REQUIRE(GM::Hydrostatic(A)() == Approx(1.0)); } SECTION("Hydrostatic - List") @@ -127,7 +126,7 @@ SECTION("Epsd - Tensor2") GM::Tensor2 A = xt::zeros({2, 2}); A(0, 1) = 1.0; A(1, 0) = 1.0; - ISCLOSE(GM::Epsd(A)(), 1.0); + REQUIRE(GM::Epsd(A)() == Approx(1.0)); } SECTION("Epsd - List") @@ -165,7 +164,7 @@ SECTION("Sigd - Tensor2") GM::Tensor2 A = xt::zeros({2, 2}); A(0, 1) = 1.0; A(1, 0) = 1.0; - ISCLOSE(GM::Sigd(A)(), 2.0); + REQUIRE(GM::Sigd(A)() == Approx(2.0)); } SECTION("Sigd - List") @@ -212,11 +211,11 @@ SECTION("Elastic - stress") mat.setStrain(Eps); auto Sig = mat.Stress(); - ISCLOSE(Sig(0, 0), K * epsm); - ISCLOSE(Sig(1, 1), K * epsm); - ISCLOSE(Sig(0, 1), G * gamma); - ISCLOSE(Sig(1, 0), G * gamma); - ISCLOSE(mat.energy(), K * std::pow(epsm, 2.0) + G * std::pow(gamma, 2.0)); + REQUIRE(Sig(0, 0) == Approx(K * epsm)); + REQUIRE(Sig(1, 1) == Approx(K * epsm)); + REQUIRE(Sig(0, 1) == Approx(G * gamma)); + REQUIRE(Sig(1, 0) == Approx(G * gamma)); + REQUIRE(mat.energy() == Approx(K * std::pow(epsm, 2.0) + G * std::pow(gamma, 2.0))); } SECTION("Cusp - stress") @@ -233,13 +232,15 @@ SECTION("Cusp - stress") mat.setStrain(Eps); auto Sig = mat.Stress(); - ISCLOSE(Sig(0, 0), K * epsm); - ISCLOSE(Sig(1, 1), K * epsm); - ISCLOSE(Sig(0, 1), 0.0); - ISCLOSE(Sig(1, 0), 0.0); - ISCLOSE(mat.epsp(), 0.02); + REQUIRE(Sig(0, 0) == Approx(K * epsm)); + REQUIRE(Sig(1, 1) == Approx(K * epsm)); + REQUIRE(Sig(0, 1) == Approx(0.0)); + REQUIRE(Sig(1, 0) == Approx(0.0)); + REQUIRE(mat.epsp() == Approx(0.02)); REQUIRE(mat.currentIndex() == 1); - ISCLOSE(mat.energy(), K * std::pow(epsm, 2.0) + G * (0.0 - std::pow(0.01, 2.0))); + REQUIRE(mat.checkYieldBoundLeft()); + REQUIRE(mat.checkYieldBoundRight()); + REQUIRE(mat.energy() == Approx(K * std::pow(epsm, 2.0) + G * (0.0 - std::pow(0.01, 2.0)))); epsm *= 2.0; gamma *= 1.9; @@ -250,13 +251,15 @@ SECTION("Cusp - stress") mat.setStrain(Eps); Sig = mat.Stress(); - ISCLOSE(Sig(0, 0), K * epsm); - ISCLOSE(Sig(1, 1), K * epsm); - ISCLOSE(Sig(0, 1), G * (gamma - 0.04)); - ISCLOSE(Sig(1, 0), G * (gamma - 0.04)); - ISCLOSE(mat.epsp(), 0.04); + REQUIRE(Sig(0, 0) == Approx(K * epsm)); + REQUIRE(Sig(1, 1) == Approx(K * epsm)); + REQUIRE(Sig(0, 1) == Approx(G * (gamma - 0.04))); + REQUIRE(Sig(1, 0) == Approx(G * (gamma - 0.04))); + REQUIRE(mat.epsp() == Approx(0.04)); REQUIRE(mat.currentIndex() == 2); - ISCLOSE(mat.energy(), K * std::pow(epsm, 2.0) + G * (std::pow(gamma - 0.04, 2.0) - std::pow(0.01, 2.0))); + REQUIRE(mat.checkYieldBoundLeft()); + REQUIRE(mat.checkYieldBoundRight()); + REQUIRE(mat.energy() == Approx(K * std::pow(epsm, 2.0) + G * (std::pow(gamma - 0.04, 2.0) - std::pow(0.01, 2.0)))); } SECTION("Smooth - stress") @@ -273,12 +276,14 @@ SECTION("Smooth - stress") mat.setStrain(Eps); auto Sig = mat.Stress(); - ISCLOSE(Sig(0, 0), K * epsm); - ISCLOSE(Sig(1, 1), K * epsm); - ISCLOSE(Sig(0, 1), 0.0); - ISCLOSE(Sig(1, 0), 0.0); - ISCLOSE(mat.epsp(), 0.02); + REQUIRE(Sig(0, 0) == Approx(K * epsm)); + REQUIRE(Sig(1, 1) == Approx(K * epsm)); + REQUIRE(Sig(0, 1) == Approx(0.0)); + REQUIRE(Sig(1, 0) == Approx(0.0)); + REQUIRE(mat.epsp() == Approx(0.02)); REQUIRE(mat.currentIndex() == 1); + REQUIRE(mat.checkYieldBoundLeft()); + REQUIRE(mat.checkYieldBoundRight()); } SECTION("Tangent (purely elastic response only)") @@ -391,6 +396,9 @@ SECTION("Array") ISCLOSE(sig(e, q, 1, 0), 0.0); ISCLOSE(epsp(e, q), fac * gamma); } + + REQUIRE(mat.checkYieldBoundLeft()); + REQUIRE(mat.checkYieldBoundRight()); } }