diff --git a/include/simcoon/Continuum_mechanics/Functions/objective_rates.hpp b/include/simcoon/Continuum_mechanics/Functions/objective_rates.hpp index dfbeaca7..92075219 100755 --- a/include/simcoon/Continuum_mechanics/Functions/objective_rates.hpp +++ b/include/simcoon/Continuum_mechanics/Functions/objective_rates.hpp @@ -152,6 +152,29 @@ void Truesdell(arma::mat &DF, arma::mat &D, arma::mat &L, const double &DTime, c */ void logarithmic_F(arma::mat &DF, arma::mat &N_1, arma::mat &N_2, arma::mat &D, arma::mat &L, const double &DTime, const arma::mat &F0, const arma::mat &F1); +/** + * @brief Computes an incremental rotation (or deformation) using the Hughes-Winget formula. + * + * This function computes an incremental rotation \f$ \Delta \mathbf{R} \f$ from a spin (or velocity gradient) tensor \f$ \mathbf{\Omega} \f$ and a time increment \f$ \Delta t \f$ + * using the Hughes-Winget midpoint formula: + * + * \f[ + * \Delta \mathbf{R} = \left( \mathbf{I} - \frac{1}{2} \Delta t \, \mathbf{\Omega} \right)^{-1} \left( \mathbf{I} + \frac{1}{2} \Delta t \, \mathbf{\Omega} \right) + * \f] + * + * @param[in] Omega 3x3 matrix representing the spin tensor \f$ \mathbf{\Omega} \f$ (or velocity gradient) + * @param[in] DTime time increment \f$ \Delta t \f$ + * @return 3x3 matrix representing the incremental rotation \f$ \Delta \mathbf{R} \f$ + * + * @details Example: + * @code + * mat Omega = randu(3,3); + * double DTime = 0.1; + * mat DR = Hughes_Winget(Omega, DTime); + * @endcode +*/ +arma::mat Hughes_Winget(const arma::mat &Omega, const double &DTime); + /** * @brief Computes the increment of rotation, the rate of deformation and the spin using the Logarithmic corotational framework using the "spin" \f$ \mathbf{\Omega}_{\textrm{log}} \f$. * @@ -378,6 +401,36 @@ arma::mat Dsigma_LieDD_2_DSDE(const arma::mat &DtauDe, const arma::mat &F); */ arma::mat DsigmaDe_JaumannDD_2_DSDE(const arma::mat &DsigmaDe, const arma::mat &F, const arma::mat &sigma); +/** + * @brief Computes the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ from the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and the approximation to logarithmic strain \f$ \mathbf{e} \f$ integrated using the Green-Naghdi spin + * + * This function takes in the tangent modulus \f$ L^t \f$ that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and the approximation to logarithmic strain \f$ \mathbf{e} \f$ integrated using the Green-Naghdi spin, + * the transformation gradient \f$ \mathbf{F} \f$ and the Kirchoff stress tensor \f$ \mathbf{\tau} \f$. + * + * It returns the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ + * + * @param[in] DtauDe (6x6 arma::mat) tangent modulus \f$ L^t \f$ that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and the approximation to logarithmic strain \f$ \mathbf{e} \f$ integrated using the Green-Naghdi spin + * @param[in] F (3x3 arma::mat) transformation gradient \f$ \mathbf{F} \f$ + * @param[in] tau (3x3 arma::mat) Kirchoff stress tensor \f$ \mathbf{\tau} \f$. + * @return (6x6 arma::mat) the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ +*/ +arma::mat DtauDe_GreenNaghdiDD_2_DSDE(const arma::mat &DtauDe, const arma::mat &F, const arma::mat &tau); + +/** + * @brief Computes the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ from the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and the approximation to logarithmic strain \f$ \mathbf{e} \f$ integrated using the Green-Naghdi spin + * + * This function takes in the tangent modulus \f$ L^t \f$ that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and the approximation to logarithmic strain \f$ \mathbf{e} \f$ integrated using the Green-Naghdi spin, + * the transformation gradient \f$ \mathbf{F} \f$ and the Cauchy stress tensor \f$ \mathbf{\sigma} \f$. + * + * It returns the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ + * + * @param[in] DsigmaDe (6x6 arma::mat) tangent modulus \f$ L^t \f$ that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and the approximation to logarithmic strain \f$ \mathbf{e} \f$ integrated using the Green-Naghdi spin + * @param[in] F (3x3 arma::mat) transformation gradient \f$ \mathbf{F} \f$ + * @param[in] sigma (3x3 arma::mat) Cauchy stress tensor \f$ \mathbf{\sigma} \f$. + * @return (6x6 arma::mat) the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ +*/ +arma::mat DsigmaDe_GreenNaghdiDD_2_DSDE(const arma::mat &DsigmaDe, const arma::mat &F, const arma::mat &sigma); + /** * @brief Computes the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and logarithmic strain \f$ \mathbf{e} \f$ (or its approximation using an integration in an rotating frame) from the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and the logarithmic strain \f$ \mathbf{e} \f$ (or its approximation using an integration in an rotating frame) * @@ -494,24 +547,75 @@ arma::mat DSDE_2_Dsigma_LieDD(const arma::mat &DSDE, const arma::mat &F); */ arma::mat DSDE_2_Dtau_JaumannDD(const arma::mat &DSDE, const arma::mat &F, const arma::mat &tau); -//arma::mat DSDE_2_Dtau_GreenNaghdiDD(const arma::mat &Lt, const arma::mat &F, const arma::mat &tau); +/** + * @brief Computes the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Green-Naghdi spin from the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ + * + * This function takes in the tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$, + * the transformation gradient \f$ \mathbf{F} \f$ and the Kirchoff stress tensor \f$ \mathbf{\tau} \f$. + * + * It returns the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Green-Naghdi spin + * + * @param[in] DSDE (6x6 arma::mat) tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ + * @param[in] F (3x3 arma::mat) transformation gradient \f$ \mathbf{F} \f$ + * @param[in] tau (3x3 arma::mat) Kirchoff stress tensor \f$ \mathbf{\tau} \f$. + * @return (6x6 arma::mat) the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Green-Naghdi spin +*/ +arma::mat DSDE_2_Dtau_GreenNaghdiDD(const arma::mat &DSDE, const arma::mat &F, const arma::mat &tau); /** - * @brief Computes the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Zaremba-Jaumann-Noll spin from the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ + * @brief Computes the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Zaremba-Jaumann-Noll spin from the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ * - * This function takes in the tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$, - * the logarithmic antisymmetric tensor-valued function \f$ \mathbf{\mathcal{B}}^{\textrm{log}} \f$, + * This function takes in the tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$, * the transformation gradient \f$ \mathbf{F} \f$ and the Cauchy stress tensor \f$ \mathbf{\sigma} \f$. - * + * * It returns the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Zaremba-Jaumann-Noll spin - * + * * @param[in] DSDE (6x6 arma::mat) tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ * @param[in] F (3x3 arma::mat) transformation gradient \f$ \mathbf{F} \f$ - * @param[in] sigma (3x3 arma::mat) Cauchy stress tensor \f$ \mathbf{\sigma} \f$. + * @param[in] sigma (3x3 arma::mat) Cauchy stress tensor \f$ \mathbf{\sigma} \f$. * @return (6x6 arma::mat) the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Zaremba-Jaumann-Noll spin */ arma::mat DSDE_2_Dsigma_JaumannDD(const arma::mat &DSDE, const arma::mat &F, const arma::mat &sigma); +/** + * @brief Computes the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Green-Naghdi spin from the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ + * + * This function takes in the tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$, + * the transformation gradient \f$ \mathbf{F} \f$ and the Cauchy stress tensor \f$ \mathbf{\sigma} \f$. + * + * It returns the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Green-Naghdi spin + * + * @param[in] DSDE (6x6 arma::mat) tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ + * @param[in] F (3x3 arma::mat) transformation gradient \f$ \mathbf{F} \f$ + * @param[in] sigma (3x3 arma::mat) Cauchy stress tensor \f$ \mathbf{\sigma} \f$. + * @return (6x6 arma::mat) the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Green-Naghdi spin +*/ +arma::mat DSDE_2_Dsigma_GreenNaghdiDD(const arma::mat &DSDE, const arma::mat &F, const arma::mat &sigma); + +/** + * @brief Computes the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the logarithmic spin from the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ + * + * Convenience function that computes get_BBBB(F) internally. + * + * @param[in] DSDE (6x6 arma::mat) tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ + * @param[in] F (3x3 arma::mat) transformation gradient \f$ \mathbf{F} \f$ + * @param[in] tau (3x3 arma::mat) Kirchoff stress tensor \f$ \mathbf{\tau} \f$. + * @return (6x6 arma::mat) the tangent modulus integrated using the logarithmic spin +*/ +arma::mat DSDE_2_Dtau_logarithmicDD(const arma::mat &DSDE, const arma::mat &F, const arma::mat &tau); + +/** + * @brief Computes the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the logarithmic spin from the tangent modulus that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$ + * + * Convenience function that computes get_BBBB(F) internally. + * + * @param[in] DSDE (6x6 arma::mat) tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ + * @param[in] F (3x3 arma::mat) transformation gradient \f$ \mathbf{F} \f$ + * @param[in] sigma (3x3 arma::mat) Cauchy stress tensor \f$ \mathbf{\sigma} \f$. + * @return (6x6 arma::mat) the tangent modulus integrated using the logarithmic spin +*/ +arma::mat DSDE_2_Dsigma_logarithmicDD(const arma::mat &DSDE, const arma::mat &F, const arma::mat &sigma); + /** * @brief Computes the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Zaremba-Jaumann-Noll spin from the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated in the natural covariant vector basis * diff --git a/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/objective_rates.cpp b/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/objective_rates.cpp index ed0276f2..8cef0880 100755 --- a/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/objective_rates.cpp +++ b/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/objective_rates.cpp @@ -50,11 +50,11 @@ py::tuple logarithmic_R(const py::array_t &F0, const py::array_t //This function computes the logarithmic strain velocity and the logarithmic spin, along with the correct rotation increment py::tuple objective_rate(const std::string& corate_name, const py::array_t &F0, const py::array_t &F1, const double &DTime, const bool &return_de, const unsigned int &n_threads) { std::map list_corate; - list_corate = { {"jaumann",0},{"green_naghdi",1},{"logarithmic",2},{"logarithmic_R",3}, {"gn",1},{"log",2},{"log_R",3}}; + list_corate = { {"jaumann",0},{"green_naghdi",1},{"logarithmic",2},{"logarithmic_R",3},{"truesdell",4},{"logarithmic_F",5}, {"gn",1},{"log",2},{"log_R",3},{"log_F",5}}; int corate = list_corate[corate_name]; - void (*corate_function)(mat &, mat &, mat &, const double &, const mat &, const mat &); - void (*corate_function_2)(mat &, mat &, mat &, mat &, mat &, const double &, const mat &, const mat &); + void (*corate_function)(mat &, mat &, mat &, const double &, const mat &, const mat &); + void (*corate_function_2)(mat &, mat &, mat &, mat &, mat &, const double &, const mat &, const mat &); switch (corate) { case 0: { @@ -72,7 +72,15 @@ py::tuple objective_rate(const std::string& corate_name, const py::array_t Delta_log_strain(const py::array_t &D, const py::arr //This function computes the logarithmic strain velocity and the logarithmic spin, along with the correct rotation increment py::array_t Lt_convert(const py::array_t &Lt, const py::array_t &F, const py::array_t &stress, const std::string &converter_key) { std::map list_Lt_convert; - list_Lt_convert = { {"Dsigma_LieDD_2_DSDE",0}, {"DsigmaDe_2_DSDE",1},{"DsigmaDe_JaumannDD_2_DSDE",2}, {"Dsigma_LieDD_Dsigma_JaumannDD",3}, {"Dsigma_LieDD_Dsigma_GreenNaghdiDD",4}, {"Dsigma_LieDD_Dsigma_logarithmicDD",5}}; + list_Lt_convert = { {"Dsigma_LieDD_2_DSDE",0}, {"DsigmaDe_2_DSDE",1},{"DsigmaDe_JaumannDD_2_DSDE",2}, {"Dsigma_LieDD_Dsigma_JaumannDD",3}, {"Dsigma_LieDD_Dsigma_GreenNaghdiDD",4}, {"Dsigma_LieDD_Dsigma_logarithmicDD",5}, {"DsigmaDe_GreenNaghdiDD_2_DSDE",6}, {"DSDE_2_Dsigma_GreenNaghdiDD",7}, {"DSDE_2_Dsigma_JaumannDD",8}, {"DSDE_2_Dsigma_LieDD",9}, {"DSDE_2_Dsigma_logarithmicDD",10}}; int select = list_Lt_convert [converter_key]; - mat (*convert_function)(const mat &, const mat &, const mat &); - mat (*convert_function2)(const mat &, const mat &); + mat (*convert_function)(const mat &, const mat &, const mat &); + mat (*convert_function2)(const mat &, const mat &); switch (select) { case 0: { - convert_function2 = &simcoon::Dsigma_LieDD_2_DSDE; + convert_function2 = &simcoon::Dsigma_LieDD_2_DSDE; break; - } + } case 1: { convert_function = &simcoon::DsigmaDe_2_DSDE; break; @@ -302,17 +330,37 @@ py::array_t Lt_convert(const py::array_t &Lt, const py::array_t< break; } case 3: { - convert_function2 = &simcoon::Dsigma_LieDD_Dsigma_JaumannDD; + convert_function2 = &simcoon::Dsigma_LieDD_Dsigma_JaumannDD; break; } case 4: { - convert_function = &simcoon::Dsigma_LieDD_Dsigma_GreenNaghdiDD; + convert_function = &simcoon::Dsigma_LieDD_Dsigma_GreenNaghdiDD; break; } case 5: { - convert_function = &simcoon::Dsigma_LieDD_Dsigma_logarithmicDD; + convert_function = &simcoon::Dsigma_LieDD_Dsigma_logarithmicDD; + break; + } + case 6: { + convert_function = &simcoon::DsigmaDe_GreenNaghdiDD_2_DSDE; break; - } + } + case 7: { + convert_function = &simcoon::DSDE_2_Dsigma_GreenNaghdiDD; + break; + } + case 8: { + convert_function = &simcoon::DSDE_2_Dsigma_JaumannDD; + break; + } + case 9: { + convert_function2 = &simcoon::DSDE_2_Dsigma_LieDD; + break; + } + case 10: { + convert_function = &simcoon::DSDE_2_Dsigma_logarithmicDD; + break; + } } if (Lt.ndim() == 2) { @@ -322,15 +370,16 @@ py::array_t Lt_convert(const py::array_t &Lt, const py::array_t< mat F_cpp = carma::arr_to_mat_view(F); mat Lt_cpp = carma::arr_to_mat_view(Lt); - vec stress_cpp = carma::arr_to_col_view(stress); + vec stress_v = carma::arr_to_col_view(stress); + mat stress_cpp = simcoon::v2t_stress(stress_v); mat Lt_converted(6,6); - switch (select) { - case 0: { + switch (select) { + case 0: case 9: { Lt_converted = convert_function2(Lt_cpp, F_cpp); break; } - case 1: case 2: { + case 1: case 2: case 6: case 7: case 8: case 10: { Lt_converted = convert_function(Lt_cpp, F_cpp, stress_cpp); break; } @@ -341,8 +390,8 @@ py::array_t Lt_convert(const py::array_t &Lt, const py::array_t< case 4: case 5: { Lt_converted = convert_function(Lt_cpp, F_cpp, stress_cpp); break; - } - } + } + } return carma::mat_to_arr(Lt_converted,false); } else if (Lt.ndim() == 3) { @@ -363,20 +412,20 @@ py::array_t Lt_convert(const py::array_t &Lt, const py::array_t< omp_set_max_active_levels(3); #pragma omp parallel for shared(Lt_converted, Lt_cpp, F_cpp) #endif -*/ +*/ for (int pt = 0; pt < nb_points; pt++) { - //vec stress_pt = stress_cpp.unsafe_col(pt); + //vec stress_pt = stress_cpp.unsafe_col(pt); stress_pt = simcoon::v2t_stress(stress_cpp.unsafe_col(pt)); - switch (select) { - case 0: { + switch (select) { + case 0: case 9: { Lt_converted.slice(pt) = convert_function2(Lt_cpp.slice(pt), F_cpp.slice(pt)); break; } - case 1: case 2: { + case 1: case 2: case 6: case 7: case 8: case 10: { Lt_converted.slice(pt) = convert_function(Lt_cpp.slice(pt), F_cpp.slice(pt), stress_pt); break; - } + } case 3: { Lt_converted.slice(pt) = convert_function2(Lt_cpp.slice(pt), stress_pt); break; @@ -384,7 +433,7 @@ py::array_t Lt_convert(const py::array_t &Lt, const py::array_t< case 4: case 5: { Lt_converted.slice(pt) = convert_function(Lt_cpp.slice(pt), F_cpp.slice(pt), stress_pt); break; - } + } } } /* #ifdef _OPENMP diff --git a/simcoon-python-builder/test/test_core/run_test.py b/simcoon-python-builder/test/test_core/run_test.py index ea022b3e..4c834c69 100755 --- a/simcoon-python-builder/test/test_core/run_test.py +++ b/simcoon-python-builder/test/test_core/run_test.py @@ -34,3 +34,118 @@ def test_isochoric_invariants(): b_vec = sim.t2v_strain(b) lambda_bar_vec = sim.isochoric_invariants(b_vec) assert np.allclose(lambda_bar, lambda_bar_vec) + + +# --------------------------------------------------------------------------- +# Shared fixture: small deformation gradient pair for objective_rate tests +# --------------------------------------------------------------------------- +@pytest.fixture(scope="module") +def deformation_pair(): + """Return (F0, F1, DTime) for a small deformation increment.""" + F0 = np.eye(3) + F1 = np.eye(3) + F1[0, 0] = 1.01 + F1[0, 1] = 0.005 + F1[1, 0] = -0.005 + F1[1, 1] = 1.02 + F1[2, 2] = 0.98 + DTime = 0.001 + return F0, F1, DTime + + +ALL_RATE_NAMES = ["jaumann", "green_naghdi", "logarithmic", "logarithmic_R", "truesdell", "logarithmic_F"] + + +def test_objective_rate_all_rates_D_consistent(deformation_pair): + """All 6 objective rates should produce the same D tensor.""" + F0, F1, DTime = deformation_pair + D_ref = None + for name in ALL_RATE_NAMES: + result = sim.objective_rate(name, F0, F1, DTime) + D = result[0] + assert D.shape == (3, 3), f"{name}: D shape {D.shape}" + # D should be symmetric + assert np.allclose(D, D.T, atol=1e-9), f"{name}: D not symmetric" + if D_ref is None: + D_ref = D + else: + assert np.allclose(D, D_ref, atol=1e-9), f"{name}: D differs from jaumann" + + +def test_objective_rate_truesdell(deformation_pair): + """Test truesdell rate: shapes and return_de.""" + F0, F1, DTime = deformation_pair + D, DR, Omega = sim.objective_rate("truesdell", F0, F1, DTime) + assert D.shape == (3, 3) + assert DR.shape == (3, 3) + assert Omega.shape == (3, 3) + # D symmetric + assert np.allclose(D, D.T, atol=1e-9) + + # return_de mode + de, D2, DR2, Omega2 = sim.objective_rate("truesdell", F0, F1, DTime, return_de=True) + assert de.shape[0] == 6 + assert np.allclose(D, D2, atol=1e-12) + + +def test_objective_rate_log_F(deformation_pair): + """Test logarithmic_F rate and its alias log_F.""" + F0, F1, DTime = deformation_pair + D1, DR1, Omega1 = sim.objective_rate("logarithmic_F", F0, F1, DTime) + D2, DR2, Omega2 = sim.objective_rate("log_F", F0, F1, DTime) + + assert D1.shape == (3, 3) + assert np.allclose(D1, D1.T, atol=1e-9) + # alias should give identical results + assert np.allclose(D1, D2, atol=1e-12) + assert np.allclose(DR1, DR2, atol=1e-12) + assert np.allclose(Omega1, Omega2, atol=1e-12) + + # return_de mode + de, D3, DR3, Omega3 = sim.objective_rate("log_F", F0, F1, DTime, return_de=True) + assert de.shape[0] == 6 + assert np.allclose(D1, D3, atol=1e-12) + + +def test_objective_rate_batch(deformation_pair): + """Test 3D batch mode (F0=2D, F1=3D) for truesdell and log_F.""" + F0, F1, DTime = deformation_pair + n_points = 4 + # Stack n_points copies of F1 into a 3×3×n cube + F1_batch = np.stack([F1] * n_points, axis=-1) + assert F1_batch.shape == (3, 3, n_points) + + for name in ["truesdell", "logarithmic_F"]: + D_batch, DR_batch, Omega_batch = sim.objective_rate(name, F0, F1_batch, DTime) + assert D_batch.shape == (3, 3, n_points), f"{name}: D_batch shape {D_batch.shape}" + assert DR_batch.shape == (3, 3, n_points) + assert Omega_batch.shape == (3, 3, n_points) + + # Single-point reference + D_single, DR_single, Omega_single = sim.objective_rate(name, F0, F1, DTime) + for i in range(n_points): + assert np.allclose(D_batch[:, :, i], D_single, atol=1e-12), f"{name}: batch D[{i}] mismatch" + + +def test_objective_rate_return_de_small_increment(): + """For very small deformation, de ≈ DTime * t2v_strain(D) for all rates.""" + F0 = np.eye(3) + eps = 1e-6 + F1 = np.eye(3) + eps * np.array([ + [1.0, 0.1, 0.0], + [0.1, -0.5, 0.0], + [0.0, 0.0, -0.5], + ]) + DTime = 1e-3 + + for name in ALL_RATE_NAMES: + result_de = sim.objective_rate(name, F0, F1, DTime, return_de=True) + de = np.asarray(result_de[0]).ravel() + D = result_de[1] + assert de.shape == (6,), f"{name}: de shape {de.shape}" + # For small increments, de ≈ DTime * t2v_strain(D) + de_approx = np.asarray(DTime * sim.t2v_strain(D)).ravel() + assert np.allclose(de, de_approx, rtol=1e-3), ( + f"{name}: de not close to DTime*t2v_strain(D), " + f"max diff = {np.max(np.abs(de - de_approx))}" + ) diff --git a/src/Continuum_mechanics/Functions/objective_rates.cpp b/src/Continuum_mechanics/Functions/objective_rates.cpp index 584ee203..ab02e901 100755 --- a/src/Continuum_mechanics/Functions/objective_rates.cpp +++ b/src/Continuum_mechanics/Functions/objective_rates.cpp @@ -251,7 +251,19 @@ void Truesdell(mat &DF, mat &D, mat &L, const double &DTime, const mat &F0, cons } catch (const std::runtime_error &e) { cerr << "Error in inv: " << e.what() << endl; throw simcoon::exception_inv("Error in inv function inside Truesdell (DF)."); - } + } +} + +mat Hughes_Winget(const mat &Omega, const double &DTime) { + mat I = eye(3,3); + mat DR; + try { + DR = (inv(I-0.5*DTime*Omega))*(I+0.5*DTime*Omega); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Hughes_Winget."); + } + return DR; } mat get_BBBB(const mat &F1) { @@ -393,7 +405,7 @@ mat DtauDe_2_DSDE(const mat &Lt, const mat &B, const mat &F, const mat &tau){ I_(i,j,k,l) = 0.5*delta_(i,k)*delta_(j,l) + 0.5*delta_(i,l)*delta_(j,k); Dtau_LieDD_(i,j,k,l) = Dtau_logarithmicDD_(i,j,k,l) + (B_(i,p,k,l)-I_(i,p,k,l))*tau_(p,j) + tau_(i,p)*(B_(j,p,k,l)-I_(j,p,k,l)); - DSDE_(L,H,M,N) = invF_(l,N)*(invF_(k,M)*(invF_(j,H)*(invF_(i,L)*Dtau_LieDD_(i,j,k,l)))); + DSDE_(L,H,M,N) = invF_(N,l)*(invF_(M,k)*(invF_(H,j)*(invF_(L,i)*Dtau_LieDD_(i,j,k,l)))); return FTensor4_mat(DSDE_); } @@ -421,7 +433,7 @@ mat Dtau_LieDD_2_DSDE(const mat &Lt, const mat &F){ Index<'M', 3> M; Index<'N', 3> N; - DSDE_(L,H,M,N) = invF_(l,N)*(invF_(k,M)*(invF_(j,H)*(invF_(i,L)*Dtau_LieDD_(i,j,k,l)))); + DSDE_(L,H,M,N) = invF_(N,l)*(invF_(M,k)*(invF_(H,j)*(invF_(L,i)*Dtau_LieDD_(i,j,k,l)))); return FTensor4_mat(DSDE_); } @@ -454,7 +466,7 @@ mat DtauDe_JaumannDD_2_DSDE(const mat &Lt, const mat &F, const mat &tau){ Index<'N', 3> N; Dtau_LieDD_(i,j,k,l) = Dtau_JaumannDD_(i,j,k,l) - 0.5*tau_(k,j)*delta_(i,l) - 0.5*tau_(l,j)*delta_(i,k) - 0.5*tau_(i,l)*delta_(j,k) - 0.5*tau_(i,k)*delta_(j,l); - DSDE_(L,H,M,N) = invF_(l,N)*(invF_(k,M)*(invF_(j,H)*(invF_(i,L)*Dtau_LieDD_(i,j,k,l)))); + DSDE_(L,H,M,N) = invF_(N,l)*(invF_(M,k)*(invF_(H,j)*(invF_(L,i)*Dtau_LieDD_(i,j,k,l)))); return FTensor4_mat(DSDE_); } @@ -468,7 +480,7 @@ mat DsigmaDe_2_DSDE(const mat &Lt, const mat &B, const mat &F, const mat &sigma) cerr << "Error in det: " << e.what() << endl; throw simcoon::exception_det("Error in det function inside DsigmaDe_2_DSDE."); } - return J*DtauDe_2_DSDE(Lt, B, F, Cauchy2Kirchoff(sigma, F, J)); + return DtauDe_2_DSDE(J*Lt, B, F, Cauchy2Kirchoff(sigma, F, J)); } //This function computes the tangent modulus that links the Piola-Kirchoff II stress S to the Green-Lagrange stress E to the tangent modulus that links the Kirchoff elastic tensor and logarithmic strain, through the log rate and the and the transformation gradient F @@ -480,35 +492,52 @@ mat DsigmaDe_2_DSDE(const mat &Lt, const mat &F, const mat &sigma){ } catch (const std::runtime_error &e) { cerr << "Error in det: " << e.what() << endl; throw simcoon::exception_det("Error in det function inside DsigmaDe_2_DSDE."); - } + } mat B = get_BBBB(F); - return J*DtauDe_2_DSDE(Lt, B, F, Cauchy2Kirchoff(sigma, F, J)); + return DtauDe_2_DSDE(J*Lt, B, F, Cauchy2Kirchoff(sigma, F, J)); } mat Dsigma_LieDD_2_DSDE(const mat &Lt, const mat &F){ - - double J; + + double J; try { J = det(F); } catch (const std::runtime_error &e) { cerr << "Error in det: " << e.what() << endl; throw simcoon::exception_det("Error in det function inside Dsigma_LieDD_2_DSDE."); - } - return J*Dtau_LieDD_2_DSDE(Lt, F); + } + return Dtau_LieDD_2_DSDE(J*Lt, F); } mat DsigmaDe_JaumannDD_2_DSDE(const mat &Lt, const mat &F, const mat &sigma){ - + double J; try { J = det(F); } catch (const std::runtime_error &e) { cerr << "Error in det: " << e.what() << endl; throw simcoon::exception_det("Error in det function inside DsigmaDe_JaumannDD_2_DSDE."); - } - return J*DtauDe_JaumannDD_2_DSDE(Lt, F, Cauchy2Kirchoff(sigma, F, J)); + } + return DtauDe_JaumannDD_2_DSDE(J*Lt, F, Cauchy2Kirchoff(sigma, F, J)); } +mat DtauDe_GreenNaghdiDD_2_DSDE(const mat &Lt, const mat &F, const mat &tau){ + + mat B = get_BBBB_GN(F); + return DtauDe_2_DSDE(Lt, B, F, tau); +} + +mat DsigmaDe_GreenNaghdiDD_2_DSDE(const mat &Lt, const mat &F, const mat &sigma){ + + double J; + try { + J = det(F); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside DsigmaDe_GreenNaghdiDD_2_DSDE."); + } + return DtauDe_GreenNaghdiDD_2_DSDE(J*Lt, F, Cauchy2Kirchoff(sigma, F, J)); +} mat DtauDe_2_DsigmaDe(const mat &Lt, const double &J) { @@ -617,17 +646,54 @@ mat DSDE_2_Dtau_JaumannDD(const mat &DSDE, const mat &F, const mat &tau) { } mat DSDE_2_Dsigma_JaumannDD(const mat &DSDE, const mat &F, const mat &sigma) { - + double J; try { J = det(F); } catch (const std::runtime_error &e) { cerr << "Error in det: " << e.what() << endl; throw simcoon::exception_det("Error in det function inside DSDE_2_Dsigma_JaumannDD."); - } + } return (1./J)*DSDE_2_Dtau_JaumannDD(DSDE, F, Cauchy2Kirchoff(sigma, F, J)); } +mat DSDE_2_Dtau_GreenNaghdiDD(const mat &DSDE, const mat &F, const mat &tau) { + + mat Dtau_LieDD = DSDE_2_Dtau_LieDD(DSDE, F); + return Dtau_LieDD_Dtau_GreenNaghdiDD(Dtau_LieDD, F, tau); +} + +mat DSDE_2_Dsigma_GreenNaghdiDD(const mat &DSDE, const mat &F, const mat &sigma) { + + double J; + try { + J = det(F); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside DSDE_2_Dsigma_GreenNaghdiDD."); + } + return (1./J)*DSDE_2_Dtau_GreenNaghdiDD(DSDE, F, Cauchy2Kirchoff(sigma, F, J)); +} + +// Standard logarithmic reverse convenience functions +mat DSDE_2_Dtau_logarithmicDD(const mat &DSDE, const mat &F, const mat &tau) { + + mat Dtau_LieDD = DSDE_2_Dtau_LieDD(DSDE, F); + return Dtau_LieDD_Dtau_logarithmicDD(Dtau_LieDD, F, tau); +} + +mat DSDE_2_Dsigma_logarithmicDD(const mat &DSDE, const mat &F, const mat &sigma) { + + double J; + try { + J = det(F); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside DSDE_2_Dsigma_logarithmicDD."); + } + return (1./J)*DSDE_2_Dtau_logarithmicDD(DSDE, F, Cauchy2Kirchoff(sigma, F, J)); +} + //This function computes the tangent modulus that links the Jaumann rate of the Kirchoff stress tau to the rate of deformation D, from the tangent modulus that links the Jaumann rate of the Kirchoff stress tau to the rate of deformation D and the Kirchoff stress tau mat Dtau_LieDD_Dtau_JaumannDD(const mat &Dtau_LieDD, const mat &tau) { diff --git a/test/Libraries/Continuum_mechanics/Tobjective_rates.cpp b/test/Libraries/Continuum_mechanics/Tobjective_rates.cpp index 97ec1b67..398911b6 100755 --- a/test/Libraries/Continuum_mechanics/Tobjective_rates.cpp +++ b/test/Libraries/Continuum_mechanics/Tobjective_rates.cpp @@ -419,3 +419,130 @@ TEST(Tobjective_rates, DSDE_conversions_LieDD) double J = det(F); EXPECT_LT(norm(Dsigma_Lie - (1./J)*Dtau_Lie, 2), 1.E-6); } + +TEST(Tobjective_rates, Hughes_Winget_identity) +{ + // Zero spin -> result should be identity + mat Omega = zeros(3,3); + double DTime = 0.001; + mat DR = Hughes_Winget(Omega, DTime); + EXPECT_LT(norm(DR - eye(3,3), 2), 1.E-12); +} + +TEST(Tobjective_rates, Hughes_Winget_orthogonal) +{ + // Antisymmetric Omega -> DR should be orthogonal (DR^T * DR = I) + mat Omega = zeros(3,3); + Omega(0,1) = 0.1; + Omega(0,2) = -0.05; + Omega(1,0) = -0.1; + Omega(1,2) = 0.08; + Omega(2,0) = 0.05; + Omega(2,1) = -0.08; + double DTime = 0.001; + + mat DR = Hughes_Winget(Omega, DTime); + EXPECT_LT(norm(DR.t() * DR - eye(3,3), 2), 1.E-6); +} + +TEST(Tobjective_rates, Hughes_Winget_consistency) +{ + // Jaumann DR should match Hughes_Winget(W, DTime) + mat F0 = eye(3,3); + mat F1 = eye(3,3); + F1(0,0) = 1.01; + F1(0,1) = 0.005; + F1(1,0) = -0.005; + F1(1,1) = 1.02; + F1(2,2) = 0.98; + double DTime = 0.001; + + mat DR_J = zeros(3,3); + mat D_J = zeros(3,3); + mat W_J = zeros(3,3); + Jaumann(DR_J, D_J, W_J, DTime, F0, F1); + + mat DR_HW = Hughes_Winget(W_J, DTime); + EXPECT_LT(norm(DR_J - DR_HW, 2), 1.E-12); +} + +TEST(Tobjective_rates, all_rates_same_D) +{ + // All 6 rate functions should produce the same D tensor + mat F0 = eye(3,3); + mat F1 = eye(3,3); + F1(0,0) = 1.01; + F1(0,1) = 0.005; + F1(1,0) = -0.005; + F1(1,1) = 1.02; + F1(2,2) = 0.98; + double DTime = 0.001; + + mat DR_j = zeros(3,3), D_j = zeros(3,3), W_j = zeros(3,3); + Jaumann(DR_j, D_j, W_j, DTime, F0, F1); + + mat DR_gn = zeros(3,3), D_gn = zeros(3,3), Omega_gn = zeros(3,3); + Green_Naghdi(DR_gn, D_gn, Omega_gn, DTime, F0, F1); + + mat DR_log = zeros(3,3), D_log = zeros(3,3), Omega_log = zeros(3,3); + logarithmic(DR_log, D_log, Omega_log, DTime, F0, F1); + + mat DR_lr = zeros(3,3), N1_lr = zeros(3,3), N2_lr = zeros(3,3), D_lr = zeros(3,3), Omega_lr = zeros(3,3); + logarithmic_R(DR_lr, N1_lr, N2_lr, D_lr, Omega_lr, DTime, F0, F1); + + mat DF_tr = zeros(3,3), D_tr = zeros(3,3), L_tr = zeros(3,3); + Truesdell(DF_tr, D_tr, L_tr, DTime, F0, F1); + + mat DF_lf = zeros(3,3), N1_lf = zeros(3,3), N2_lf = zeros(3,3), D_lf = zeros(3,3), L_lf = zeros(3,3); + logarithmic_F(DF_lf, N1_lf, N2_lf, D_lf, L_lf, DTime, F0, F1); + + EXPECT_LT(norm(D_j - D_gn, 2), 1.E-9); + EXPECT_LT(norm(D_j - D_log, 2), 1.E-9); + EXPECT_LT(norm(D_j - D_lr, 2), 1.E-9); + EXPECT_LT(norm(D_j - D_tr, 2), 1.E-9); + EXPECT_LT(norm(D_j - D_lf, 2), 1.E-9); +} + +TEST(Tobjective_rates, Truesdell_logF_same_DF) +{ + // Truesdell and logarithmic_F should produce the same DF (both use Hughes-Winget on L) + mat F0 = eye(3,3); + mat F1 = eye(3,3); + F1(0,0) = 1.01; + F1(0,1) = 0.005; + F1(1,0) = -0.005; + F1(1,1) = 1.02; + F1(2,2) = 0.98; + double DTime = 0.001; + + mat DF_tr = zeros(3,3), D_tr = zeros(3,3), L_tr = zeros(3,3); + Truesdell(DF_tr, D_tr, L_tr, DTime, F0, F1); + + mat DF_lf = zeros(3,3), N1_lf = zeros(3,3), N2_lf = zeros(3,3), D_lf = zeros(3,3), L_lf = zeros(3,3); + logarithmic_F(DF_lf, N1_lf, N2_lf, D_lf, L_lf, DTime, F0, F1); + + EXPECT_LT(norm(DF_tr - DF_lf, 2), 1.E-12); + EXPECT_LT(norm(L_tr - L_lf, 2), 1.E-12); +} + +TEST(Tobjective_rates, logR_logF_same_N1) +{ + // logarithmic_R and logarithmic_F should produce the same N_1 (same f_log spectral function) + mat F0 = eye(3,3); + mat F1 = eye(3,3); + F1(0,0) = 1.01; + F1(0,1) = 0.005; + F1(1,0) = -0.005; + F1(1,1) = 1.02; + F1(2,2) = 0.98; + double DTime = 0.001; + + mat DR_lr = zeros(3,3), N1_lr = zeros(3,3), N2_lr = zeros(3,3), D_lr = zeros(3,3), Omega_lr = zeros(3,3); + logarithmic_R(DR_lr, N1_lr, N2_lr, D_lr, Omega_lr, DTime, F0, F1); + + mat DF_lf = zeros(3,3), N1_lf = zeros(3,3), N2_lf = zeros(3,3), D_lf = zeros(3,3), L_lf = zeros(3,3); + logarithmic_F(DF_lf, N1_lf, N2_lf, D_lf, L_lf, DTime, F0, F1); + + EXPECT_LT(norm(N1_lr - N1_lf, 2), 1.E-12); + EXPECT_LT(norm(N2_lr - N2_lf, 2), 1.E-12); +}