diff --git a/CMakeLists.txt b/CMakeLists.txt index d229db02e..0c351c6be 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -357,11 +357,17 @@ if(SIMCOON_BUILD_TESTS) enable_language(Fortran) add_library(umat_plugin_aba SHARED testBin/Umats/UMABA/external/UMAT_ABAQUS_ELASTIC.for testBin/Umats/UMABA/external/umat_plugin_aba.cpp) - set_target_properties(umat_plugin_aba PROPERTIES - LIBRARY_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/testBin/Umats/UMABA/external - ) +if(NOT MSVC) + set_source_files_properties(testBin/Umats/UMABA/external/UMAT_ABAQUS_ELASTIC.f PROPERTIES Fortran_PREPROCESS OFF) + add_library(umat_plugin_aba SHARED testBin/Umats/UMABA/external/UMAT_ABAQUS_ELASTIC.f testBin/Umats/UMABA/external/umat_plugin_aba.cpp) + set_target_properties(umat_plugin_aba PROPERTIES + PREFIX "" + SUFFIX "" + LIBRARY_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/testBin/Umats/UMABA/external + ) - target_link_libraries(umat_plugin_aba PRIVATE simcoon) + target_link_libraries(umat_plugin_aba PRIVATE simcoon ${ARMADILLO_LIBRARIES}) +endif() endif() add_library(umat_plugin_ext SHARED diff --git a/environment_arm64.yml b/environment_arm64.yml index 14f3eb304..3605c4965 100644 --- a/environment_arm64.yml +++ b/environment_arm64.yml @@ -6,7 +6,7 @@ channels: dependencies: - cxx-compiler - - fortran-compiler + - gfortran_osx-arm64 - blas - lapack - numpy>=1.26 diff --git a/external/UMAT_ABAQUS_ELASTIC.for b/external/UMAT_ABAQUS_ELASTIC.f similarity index 100% rename from external/UMAT_ABAQUS_ELASTIC.for rename to external/UMAT_ABAQUS_ELASTIC.f diff --git a/include/simcoon/Continuum_mechanics/Functions/contimech.hpp b/include/simcoon/Continuum_mechanics/Functions/contimech.hpp index 5536c673e..8e58872c4 100755 --- a/include/simcoon/Continuum_mechanics/Functions/contimech.hpp +++ b/include/simcoon/Continuum_mechanics/Functions/contimech.hpp @@ -492,7 +492,7 @@ arma::mat sym_dyadic_operator(const arma::mat &a, const arma::mat &b); double f_z = 0.; for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { - if ((i!=j)&&(fabs(bi(i)-bi(j))>sim_iota)) { + if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) { f_z = (1.+(bi(i)/bi(j)))/(1.-(bi(i)/bi(j)))+2./log(bi(i)/bi(j)); BBBB = BBBB + f_z*B_klmn(Bi.col(i),Bi.col(j)); } diff --git a/include/simcoon/Continuum_mechanics/Functions/recovery_props.hpp b/include/simcoon/Continuum_mechanics/Functions/recovery_props.hpp index 9cff4e668..9c1c73950 100755 --- a/include/simcoon/Continuum_mechanics/Functions/recovery_props.hpp +++ b/include/simcoon/Continuum_mechanics/Functions/recovery_props.hpp @@ -42,7 +42,7 @@ namespace simcoon{ * @param[out] axis : The axis of symmetry (if applicable) * @param[out] props : The material properties vector * @param[out] maj_sym : The major symmetry condition (L_ij = L_ji ?). - * @param[in] tol : The tolerance utilized to check the symetries. If less than the global sim_limit (1.E-8), sim_limit is utilized. Default is 0. (so sim_limit is utilized by default) + * @param[in] tol : The tolerance utilized to check the symetries. If less than the global simcoon::limit (1.E-8), simcoon::limit is utilized. Default is 0. (so simcoon::limit is utilized by default) * * Material Symmetries considered: * diff --git a/include/simcoon/Continuum_mechanics/Material/ODF.hpp b/include/simcoon/Continuum_mechanics/Material/ODF.hpp index 0f71b262b..4bbfa7f55 100755 --- a/include/simcoon/Continuum_mechanics/Material/ODF.hpp +++ b/include/simcoon/Continuum_mechanics/Material/ODF.hpp @@ -62,7 +62,7 @@ class ODF arma::vec limits; //minimal and maximal angles of orientation ODF(); //default constructor - ODF(const int &, const bool & = false, const double & = 0., const double & = sim_pi); //Partial constructor + ODF(const int &, const bool & = false, const double & = 0., const double & = simcoon::pi); //Partial constructor ODF(const int &, const int &, const int &, const std::vector &, const double &, const arma::vec &, const arma::vec &, const bool & = false); //Full constructor ODF(const ODF&); //Copy constructor diff --git a/include/simcoon/exception.hpp b/include/simcoon/exception.hpp new file mode 100755 index 000000000..0d91f1081 --- /dev/null +++ b/include/simcoon/exception.hpp @@ -0,0 +1,71 @@ +/* This file is part of simcoon. + + simcoon is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + simcoon is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with simcoon. If not, see . + + */ + +///@file exception.hpp +///@brief exceptions of simcoon +///@version 1.0 + +#pragma once + +#include +#include + +namespace simcoon { + +class exception_eig_sym : public std::runtime_error { +public: + explicit exception_eig_sym(const std::string& msg) + : std::runtime_error(msg) {} +}; + +class exception_det : public std::runtime_error { + public: + explicit exception_det(const std::string& msg) + : std::runtime_error(msg) {} + }; + +class exception_inv : public std::runtime_error { + public: + explicit exception_inv(const std::string& msg) + : std::runtime_error(msg) {} + }; + +class exception_sqrtmat_sympd : public std::runtime_error { + public: + explicit exception_sqrtmat_sympd(const std::string& msg) + : std::runtime_error(msg) {} + }; + +class exception_logmat_sympd : public std::runtime_error { + public: + explicit exception_logmat_sympd(const std::string& msg) + : std::runtime_error(msg) {} + }; + +class exception_expmat_sym : public std::runtime_error { + public: + explicit exception_expmat_sym(const std::string& msg) + : std::runtime_error(msg) {} + }; + +class exception_powmat : public std::runtime_error { + public: + explicit exception_powmat(const std::string& msg) + : std::runtime_error(msg) {} + }; + +} //namespace simcoon diff --git a/include/simcoon/parameter.hpp b/include/simcoon/parameter.hpp index 5280638be..323bf76c9 100755 --- a/include/simcoon/parameter.hpp +++ b/include/simcoon/parameter.hpp @@ -19,6 +19,7 @@ ///@brief parameters of simcoon ///@version 1.0 +#pragma once #if defined(__cpp_lib_math_constants) #include #define HAS_STD_NUMBERS @@ -26,66 +27,29 @@ #define UNUSED(x) [&x]{}() -#ifndef version_full -#define version_full 1 -#endif - namespace simcoon{ -#ifndef sim_pi #ifdef HAS_STD_NUMBERS - #define sim_pi std::numbers::pi +constexpr double pi = std::numbers::pi; #else - #define sim_pi 3.14159265358979323846 +constexpr double pi = 3.14159265358979323846; #endif -#endif -#ifndef axis_psi -#define axis_psi 3 -#endif +constexpr int axis_psi = 3; +constexpr int axis_theta = 1; +constexpr int axis_phi = 3; -#ifndef axis_theta -#define axis_theta 1 -#endif +constexpr double limit = 1.E-9; +constexpr double iota = 1.E-12; +constexpr int miniter_umat = 10; -#ifndef axis_phi -#define axis_phi 3 -#endif +constexpr int maxiter_umat = 100; +constexpr double precision_umat = 1E-9; -#ifndef sim_limit -#define sim_limit 1.E-9 -#endif +constexpr double div_tnew_dt_umat = 0.2; +constexpr double mul_tnew_dt_umat = 2.0; +constexpr int maxiter_micro = 100; -#ifndef sim_iota -#define sim_iota 1.E-12 -#endif - -#ifndef miniter_umat -#define miniter_umat 10 -#endif - -#ifndef maxiter_umat -#define maxiter_umat 100 -#endif - -#ifndef precision_umat -#define precision_umat 1E-9 -#endif - -#ifndef div_tnew_dt_umat -#define div_tnew_dt_umat 0.2 -#endif - -#ifndef mul_tnew_dt_umat -#define mul_tnew_dt_umat 2 -#endif - -#ifndef maxiter_micro -#define maxiter_micro 100 -#endif - -#ifndef precision_micro -#define precision_micro 1E-6 -#endif +constexpr double precision_micro = 1E-6; } //namespace simcoon diff --git a/simcoon-python-builder/CMakeLists.txt b/simcoon-python-builder/CMakeLists.txt index 78302a69a..557975a3c 100755 --- a/simcoon-python-builder/CMakeLists.txt +++ b/simcoon-python-builder/CMakeLists.txt @@ -13,7 +13,147 @@ else() message(STATUS "Building Python bindings (standalone)") endif() -# All dependencies (Python, pybind11, carma, Armadillo, simcoon) are already found by parent +#set(CMAKE_VERBOSE_MAKEFILE ON) + +# The version number +set(VERSION_MAJOR 1) +set(VERSION_MINOR 9) +set(VERSION_PATCH 4) + +if(VERSION_PATCH MATCHES "0") + set(VERSION_NUMBER "${VERSION_MAJOR}.${VERSION_MINOR}") +else() + set(VERSION_NUMBER "${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}") +endif() + +set(CMAKE_DISABLE_PRECOMPILE_HEADERS ON) + +##### RPATH HANDLING #### +#set(CMAKE_MACOSX_RPATH 1) +#set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") +#set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) +################ + + +# put our local cmake find scripts at the beginning of the cmake +# module search path +set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH}) + +set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) +set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) + +if (MSVC) #if windows using visual c++ + + # Specify the output directory for each configuration (Debug, Release, etc.) + foreach(CONFIG_TYPE ${CMAKE_CONFIGURATION_TYPES}) + string(TOUPPER ${CONFIG_TYPE} CONFIG) + set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY_${CONFIG} "${CMAKE_BINARY_DIR}/lib") + set(CMAKE_LIBRARY_OUTPUT_DIRECTORY_${CONFIG} "${CMAKE_BINARY_DIR}/lib") + set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_${CONFIG} "${CMAKE_BINARY_DIR}/bin") + endforeach() +endif() + +#Set executable files and library files +set(EXECUTABLE_OUTPUT_PATH bin/${CMAKE_BUILD_TYPE}) +set(LIBRARY_OUTPUT_PATH lib/${CMAKE_BUILD_TYPE}) + +find_package(Python3 COMPONENTS Interpreter Development NumPy) + +#Inclusion of Armadillo +if (MSVC) + find_package(pybind11 REQUIRED) + find_package(carma CONFIG REQUIRED) + find_package(blas REQUIRED) + find_package(lapack REQUIRED) + set(ARMADILLO_LIBRARIES ${ARMADILLO_LIBRARIES} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES}) +else() + find_package(pybind11 REQUIRED) + find_package(Armadillo 12.6 REQUIRED) + find_package(carma CONFIG REQUIRED) +endif() + +if (NOT DEFINED USE_OPENMP) + set(USE_OPENMP True) +endif() + +if (USE_OPENMP) + find_package(OpenMP) +endif() +find_package(simcoon REQUIRED) + +#${PYTHON_VERSION_MAJOR}${PYTHON_VERSION_MINOR} + +include_directories(SYSTEM ${ARMADILLO_INCLUDE_DIRS}) +include_directories(SYSTEM ${PYTHON_INCLUDE_DIRS}) +include_directories(SYSTEM ${SIMCOON_INCLUDE_DIRS}) +include_directories(SYSTEM ${CARMA_INCLUDE_DIR}) + +message(STATUS "PYTHON_LIBRARIES = ${PYTHON_LIBRARIES}") +message(STATUS "PYTHON_EXECUTABLE = ${PYTHON_EXECUTABLE}") +message(STATUS "PYTHON_INCLUDE_DIRS = ${PYTHON_INCLUDE_DIRS}") +message(STATUS "SIMCOON_INCLUDE_DIRS = ${SIMCOON_INCLUDE_DIRS}") +message(STATUS "SIMCOON_LIBRARIES = ${SIMCOON_LIBRARIES}") +message(STATUS "ARMADILLO_INCLUDE_DIRS = ${ARMADILLO_INCLUDE_DIRS}") +message(STATUS "ARMADILLO_LIBRARIES = ${ARMADILLO_LIBRARIES}") +message(STATUS "carma_INCLUDE_DIR = ${carma_INCLUDE_DIR}") + +# Build type +if(NOT CMAKE_BUILD_TYPE) # Debug by default + set(CMAKE_BUILD_TYPE Debug CACHE STRING + "Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel" + FORCE) +endif() + +# Set the CFLAGS and CXXFLAGS depending on the options the user specified. +#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra") +#set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wall -Wextra -Wno-unused-parameter") + + +if (MSVC) + set(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS ON) + if (USE_OPENMP) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /std:c++latest /Drestrict= /openmp:llvm /Y-") #\Y- to disable precompile header (don't work for an unknown reason) + # set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /std:c++latest /Drestrict= /openmp /Y-") #\Y- to disable precompile header (don't work for an unknown reason) + else() + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /std:c++latest /Drestrict=") + endif() + set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS}") #do nothing but kept anyway if required to add some options + set(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS}") #do nothing but kept anyway if required to add some options + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") #do nothing but kept anyway if required to add some options +elseif (UNIX AND NOT APPLE) + if (DEBUG) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++20 -O0 -Drestrict=") + set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} -std=c++20") + set(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -std=c++20") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -std=c++20") + else() + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++20 -O3 -Drestrict=") + set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} -std=c++20") + set(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -std=c++20") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -std=c++20") + endif() +else() + if (DEBUG) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++20 -stdlib=libc++ -O0 -Drestrict=") + set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} -std=c++20 -stdlib=libc++") + set(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -std=c++20 -stdlib=libc++") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -std=c++20 -stdlib=libc++") + else() + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++20 -stdlib=libc++ -O3 -Drestrict=") + set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} -std=c++20 -stdlib=libc++") + set(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -std=c++20 -stdlib=libc++") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -std=c++20 -stdlib=libc++") + endif() +endif() + +#Inclusion of public headers +include_directories(include) + +#Command file to get all the files in the src/ and include/ +file(GLOB_RECURSE source_files_simmit src/python_wrappers/* include/simcoon/python_wrappers/*) +file(GLOB_RECURSE test_files_carma test/Carma/*.cpp test/Carma/*.hpp) +file(GLOB_RECURSE test_files_exception test/exception_test/*.cpp test/exception_test/*.hpp) ################################################################################ # SIMCOON PYBIND11 MODULE (simmit) diff --git a/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/hyperelastic.cpp b/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/hyperelastic.cpp index dc69b0bae..7a1a24e1a 100755 --- a/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/hyperelastic.cpp +++ b/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/hyperelastic.cpp @@ -4,8 +4,7 @@ #include #include - - +#include #include #include #include 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 862b987c6..ed0276f20 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 @@ -10,6 +10,7 @@ #include #endif +#include #include #include #include @@ -84,7 +85,7 @@ py::tuple objective_rate(const std::string& corate_name, const py::array_t sim_limit) { + if (Time > simcoon::limit) { start = false; } diff --git a/simcoon-python-builder/src/python_wrappers/boostpython_smartplus_wrappers.cpp b/simcoon-python-builder/src/python_wrappers/boostpython_smartplus_wrappers.cpp index 467e96609..053c91ff3 100755 --- a/simcoon-python-builder/src/python_wrappers/boostpython_smartplus_wrappers.cpp +++ b/simcoon-python-builder/src/python_wrappers/boostpython_smartplus_wrappers.cpp @@ -3,6 +3,8 @@ #include #include +#include + #include #include #include @@ -54,6 +56,18 @@ PYBIND11_MODULE(simmit, m) { m.doc() = "pybind11 example plugin"; // optional module docstring + // Create a Python-visible base exception for all your custom errors + py::object SimcoonError = py::register_exception(m, "SimcoonError", PyExc_RuntimeError); + + // Register the exception translator + py::register_exception(m, "CppExceptionEigSym", SimcoonError.ptr()); + py::register_exception(m, "CppExceptionDet", SimcoonError.ptr()); + py::register_exception(m, "CppExceptionInv", SimcoonError.ptr()); + py::register_exception(m, "CppExceptionSqrtMatSym", SimcoonError.ptr()); + py::register_exception(m, "CppExceptionLogMatSym", SimcoonError.ptr()); + py::register_exception(m, "CppExceptionExpMatSym", SimcoonError.ptr()); + py::register_exception(m, "CppExceptionPowMat", SimcoonError.ptr()); + // Register the from-python converters for constitutive.hpp m.def("Ireal", &Ireal, "copy"_a = true, simcoon_docs::Ireal); m.def("Ivol", &Ivol, "copy"_a = true, simcoon_docs::Ivol); diff --git a/simcoon-python-builder/test/exception_test/Teig_sym.cpp b/simcoon-python-builder/test/exception_test/Teig_sym.cpp new file mode 100644 index 000000000..207cc4c73 --- /dev/null +++ b/simcoon-python-builder/test/exception_test/Teig_sym.cpp @@ -0,0 +1,88 @@ +/* This file is part of simcoon. + + simcoon is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + simcoon is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with simcoon. If not, see . + + */ + +///@file Teig_sym.cpp +///@brief Test for the exceptions raised with simcoon +///@version 1.0 + +#include +#include +#include +#include +#include +#include + +using namespace std; +namespace py = pybind11; +namespace testexception { + +py::array_t test_eig_sym_val_affect(const py::array_t &matrix_fail_eigen, const bool ©) { + + //Try eigen decomposition : fail since decomposition is not unique + arma::mat arma_fail_eigen = carma::arr_to_mat(matrix_fail_eigen); + + arma::vec eigval; + + try { + eigval = arma::eig_sym(arma_fail_eigen); + } catch (const std::runtime_error &e) { + cerr << "Error in eig_sym: " << e.what() << endl; + throw simcoon::exception_eig_sym("Error in eig_sym function inside test_eig_sym_val_affect."); + } +} + +py::array_t test_eig_sym_val_modify(const py::array_t &matrix_fail_eigen, const bool ©) { + + //Try eigen decomposition : fail since decomposition is not unique + arma::mat arma_fail_eigen = carma::arr_to_mat(matrix_fail_eigen); + arma::vec eigval; + + bool success = arma::eig_sym(eigval, arma_fail_eigen); + + if (success) { + eigval.print("Eigenvalues:"); + return carma::col_to_arr(eigval, copy); + } else { + cout << "Eigen decomposition failed." << endl; + throw simcoon::exception_eig_sym("Error in eig_sym function inside test_eig_sym_val_modify."); + } +} + +py::tuple test_eig_sym_val_vec(const py::array_t &matrix_fail_eigen, const bool ©) { + + //Try eigen decomposition : fail since decomposition is not unique + arma::mat arma_fail_eigen = carma::arr_to_mat(matrix_fail_eigen); + arma::vec eigval; + arma::mat eigvec; + + bool success = arma::eig_sym(eigval, eigvec, arma_fail_eigen); + + if (success) { + eigval.print("Eigenvalues:"); + eigvec.print("Eigenvectors:"); + return py::make_tuple(carma::col_to_arr(eigval, copy), carma::mat_to_arr(eigvec, copy)); + } else { + cout << "Eigen decomposition failed." << endl; + throw simcoon::exception_eig_sym("Error in eig_sym function inside test_eig_sym_val_vec."); + } +} + + + + + +} // namespace testexception diff --git a/simcoon-python-builder/test/exception_test/Teig_sym.hpp b/simcoon-python-builder/test/exception_test/Teig_sym.hpp new file mode 100644 index 000000000..c7d3c5d93 --- /dev/null +++ b/simcoon-python-builder/test/exception_test/Teig_sym.hpp @@ -0,0 +1,38 @@ +/* This file is part of simcoon. + + simcoon is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + simcoon is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with arma2numpy. If not, see . + + */ + +///@file Teig_sym.hpp +///@brief Test for the exceptions raised with simcoon +///@version 1.0 + +#pragma once +#include +#include +#include +#include + +namespace testexception { + +pybind11::array_t test_eig_sym_val_affect(const py::array_t &matrix_fail_eigen, const bool ©=true); + +pybind11::array_t test_eig_sym_val_modify(const py::array_t &matrix_fail_eigen, const bool ©=true); + +pybind11::tuple test_eig_sym_val_vec(const py::array_t &matrix_fail_eigen, const bool ©=true); + +//py::tuple test_eig_sym(py::array_t& arr, bool copy); + +} // namespace testexception \ No newline at end of file diff --git a/simcoon-python-builder/test/exception_test/bindings.cpp b/simcoon-python-builder/test/exception_test/bindings.cpp new file mode 100755 index 000000000..d53f02628 --- /dev/null +++ b/simcoon-python-builder/test/exception_test/bindings.cpp @@ -0,0 +1,25 @@ +#include +// include numpy header for usage of array_t +#include + +#include +#include +#include +#include "Teig_sym.hpp" + +#include + +namespace py = pybind11; +using namespace testexception; + +using namespace pybind11::literals; + +PYBIND11_MODULE(test_exception, m) { + + //Register the test_eig_sym_val_vec + m.def("test_eig_sym_val_modify", &test_eig_sym_val_modify, "matrix_fail_eigen"_a, "copy"_a=true, "This function is used to test the eig_sym exception (modify eigval)"); + m.def("test_eig_sym_val_affect", &test_eig_sym_val_affect, "matrix_fail_eigen"_a, "copy"_a=true, "This function is used to test the eig_sym exception (affectation)"); + m.def("test_eig_sym_val_vec", &test_eig_sym_val_vec, "matrix_fail_eigen"_a, "copy"_a=true, "This function is used to test the eig_sym exception (modify eigval and eigvec)"); + + +} diff --git a/simcoon-python-builder/test/exception_test/run_test_exception.py b/simcoon-python-builder/test/exception_test/run_test_exception.py new file mode 100644 index 000000000..834ce5bf3 --- /dev/null +++ b/simcoon-python-builder/test/exception_test/run_test_exception.py @@ -0,0 +1,48 @@ +"""Test numpy array to matrix conversion function.""" + +import numpy as np +import numpy.typing as npt +import test_exception as exception +from simcoon import simmit as sim +import pytest + + +@pytest.fixture(scope="session") +def Mat_fail_eigen_sym() -> npt.NDArray[np.float64]: + mat_fail_eigen_sym = np.array( + [ + [1.0e1222, 2, 3], + [9, 1, 4], # 9 ≠ 2 → asymmetric in lower triangle + [0, 0, 1], + ] + ) + return mat_fail_eigen_sym + + +def test_eig_sym(Mat_fail_eigen_sym): + """Test eigen_sym exception.""" + + try: + eigval = exception.test_eig_sym_val_affect(Mat_fail_eigen_sym, True) + except sim.SimcoonError as e: + print("Caught C++ exception in test_eig_sym_val:", e) + + try: + eigval = exception.test_eig_sym_val_affect(Mat_fail_eigen_sym, True) + except sim.SimcoonError as e: + print("Caught C++ exception in test_eig_sym_val:", e) + + try: + eigval = exception.test_eig_sym_val_affect(Mat_fail_eigen_sym, True) + except RuntimeError as e: + print("Caught C++ std::runtime_error exception in test_eig_sym_val:", e) + + try: + eigval = exception.test_eig_sym_val_modify(Mat_fail_eigen_sym, True) + except sim.SimcoonError as e: + print("Caught C++ exception in test_eig_sym_val_vec:", e) + + try: + (eigval, eigvec) = exception.test_eig_sym_val_vec(Mat_fail_eigen_sym, True) + except sim.SimcoonError as e: + print("Caught C++ exception in test_eig_sym_val_vec:", e) diff --git a/simcoon-python-builder/test/exception_test/test_exception.so b/simcoon-python-builder/test/exception_test/test_exception.so new file mode 100755 index 000000000..1edf4ceb5 Binary files /dev/null and b/simcoon-python-builder/test/exception_test/test_exception.so differ diff --git a/software/ODF.cpp b/software/ODF.cpp index 1e13f4af9..1569c1993 100755 --- a/software/ODF.cpp +++ b/software/ODF.cpp @@ -107,7 +107,7 @@ int main() { string npeaksfile = "Npeaks" + to_string(int(num_file_npeaks)) + ".dat"; read_peak(odf_rve, path_data, npeaksfile); - vec x = linspace(angle_min, angle_max-sim_iota, 90); + vec x = linspace(angle_min, angle_max-simcoon::iota, 90); cout << "x = " << x.t() << endl; vec y = get_densities_ODF(x, path_data, "Npeaks0.dat", false); diff --git a/software/PDF.cpp b/software/PDF.cpp index e2f0ee489..13c86eb4f 100755 --- a/software/PDF.cpp +++ b/software/PDF.cpp @@ -103,7 +103,7 @@ int main() { string npeaksfile = "Npeaks" + to_string(int(num_file_npeaks)) + ".dat"; read_peak(pdf_rve, path_data, npeaksfile); - vec x = linspace(parameter_min, parameter_max-sim_iota, 100 ); + vec x = linspace(parameter_min, parameter_max-simcoon::iota, 100 ); cout << "x = " << x.t() << endl; vec y = get_densities_PDF(x, path_data, npeaksfile); diff --git a/src/Continuum_mechanics/Functions/constitutive.cpp b/src/Continuum_mechanics/Functions/constitutive.cpp index 6c9b7f2be..be3f73743 100755 --- a/src/Continuum_mechanics/Functions/constitutive.cpp +++ b/src/Continuum_mechanics/Functions/constitutive.cpp @@ -24,6 +24,7 @@ along with simcoon. If not, see . #include #include #include +#include #include using namespace std; @@ -384,7 +385,12 @@ mat M_ortho(const double &C11, const double &C12, const double &C13, const doubl L(4,4) = C55; L(5,5) = C66; - M = inv(L); + try { + M = inv(L); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside M_ortho."); + } } else if (conv == "EnuG") { diff --git a/src/Continuum_mechanics/Functions/contimech.cpp b/src/Continuum_mechanics/Functions/contimech.cpp index 2c3f1a626..91bc584ea 100755 --- a/src/Continuum_mechanics/Functions/contimech.cpp +++ b/src/Continuum_mechanics/Functions/contimech.cpp @@ -257,9 +257,9 @@ double Macaulay_n(const double &d) { //This function returns the value if it's negative, zero if it's positive (Macaulay brackets <>-) double sign(const double &d) { - if ((d < sim_iota)&&(fabs(d) > sim_iota)) + if ((d < simcoon::iota)&&(fabs(d) > simcoon::iota)) return -1.; - else if(d > sim_iota) + else if(d > simcoon::iota) return 1.; else return 0.; @@ -355,28 +355,16 @@ mat p_ikjl(const vec &a) { mat auto_sym_dyadic(const mat &A) { //T vec A_v = t2v_sym(A); - mat C = zeros(6,6); - for (unsigned int i=0; i<6; i++) { - for (unsigned int j=0; j<6; j++) { - C(i,j) = A_v(i)*A_v(j); - } - } - return C; + return A_v * A_v.t(); } mat sym_dyadic(const mat &A, const mat &B) { vec A_v = t2v_sym(A); vec B_v = t2v_sym(B); - mat C = zeros(6,6); - - for (unsigned int i=0; i<6; i++) { - for (unsigned int j=0; j<6; j++) { - C(i,j) = A_v(i)*B_v(j); - } - } - return C; + + return A_v * B_v.t(); } mat auto_dyadic(const mat &A) { diff --git a/src/Continuum_mechanics/Functions/criteria.cpp b/src/Continuum_mechanics/Functions/criteria.cpp index e4048587d..ba42f535f 100755 --- a/src/Continuum_mechanics/Functions/criteria.cpp +++ b/src/Continuum_mechanics/Functions/criteria.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include diff --git a/src/Continuum_mechanics/Functions/derivatives.cpp b/src/Continuum_mechanics/Functions/derivatives.cpp index 494536062..465646a3d 100755 --- a/src/Continuum_mechanics/Functions/derivatives.cpp +++ b/src/Continuum_mechanics/Functions/derivatives.cpp @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -58,13 +59,26 @@ mat dtrSdS(const mat &S) { //This function returns the derivative of the determinant of a tensor mat ddetSdS(const mat &S) { - return det(S)*(inv(S)).t(); + + try { + return det(S)*(inv(S)).t(); + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv (combined expression), error inv is thrown: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside ddetSdS."); + } } mat dinvSdSsym(const mat &S) { - mat invS = inv(S); - + mat invS; + + try { + invS = inv(S); + } catch (const std::runtime_error &e) { + cerr << "Error in inv : " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside dinvSdSsym."); + } + Tensor2 invS_ = mat_FTensor2(invS); Tensor4 dinvSdSsym_; diff --git a/src/Continuum_mechanics/Functions/hyperelastic.cpp b/src/Continuum_mechanics/Functions/hyperelastic.cpp index b4f4be6c1..5558ea271 100755 --- a/src/Continuum_mechanics/Functions/hyperelastic.cpp +++ b/src/Continuum_mechanics/Functions/hyperelastic.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -33,13 +34,24 @@ namespace simcoon{ vec isochoric_invariants(const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = sqrt(det(b)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(b)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside isochoric_invariants."); + } } mat b_bar = pow(J,-2./3.)*b; vec I = zeros(3); + I(0) = trace(b_bar); - I(1) = 0.5*(pow(trace(b_bar),2.)-trace(powmat(b_bar,2))); + try { + I(1) = 0.5*(pow(trace(b_bar),2.)-trace(powmat(b_bar,2))); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_powmat("Error in powmat function inside isochoric_invariants."); + } I(2) = 1.; return I; } @@ -47,7 +59,7 @@ vec isochoric_invariants(const mat &b, const double &mJ) { vec isochoric_invariants(const vec &lambda, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { + if (fabs(mJ) < simcoon::iota) { J = prod(lambda); } vec lambda_bar = pow(J,-1./3.)*lambda; @@ -60,10 +72,21 @@ vec isochoric_invariants(const vec &lambda, const double &mJ) { vec isochoric_pstretch_from_V(const mat &V, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = det(V); + if (fabs(mJ) < simcoon::iota) { + try { + J = det(V); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside isochoric_pstretch_from_V."); + } + } + vec lambda; + try { + lambda = eig_sym(V); + } catch (const std::runtime_error &e) { + cerr << "Error in eig_sym: " << e.what() << endl; + throw simcoon::exception_eig_sym("Failed to compute eigenvalues in isochoric_pstretch_from_V."); } - vec lambda = eig_sym(V); vec lambda_bar = pow(J,-1./3.)*lambda; return lambda_bar; } @@ -71,10 +94,21 @@ vec isochoric_pstretch_from_V(const mat &V, const double &mJ) { vec isochoric_pstretch_from_b(const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = sqrt(det(b)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(b)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside isochoric_pstretch_from_b."); + } + } + vec lambda; + try { + lambda = eig_sym(b); + } catch (const std::runtime_error &e) { + cerr << "Error in eig_sym: " << e.what() << endl; + throw simcoon::exception_eig_sym("Failed to compute eigenvalues in isochoric_pstretch_from_b."); } - vec lambda = eig_sym(b); lambda.transform( [](double val) { return (sqrt(val)); } ); vec lambda_bar = pow(J,-1./3.)*lambda; return lambda_bar; @@ -90,33 +124,49 @@ vec isochoric_pstretch(const mat &input, const string &input_tensor, const doubl throw std::invalid_argument("Invalid input string to describe the input vector: it should be *b* for left Cauchy-Green tensor or *v* or *V* for Eulerian stretch tensor"); } -void pstretch(vec &lambda, mat &n_pvector, const mat &input, const string &input_tensor, const double &mJ) { +void pstretch(vec &lambda, mat &n_pvectors, const mat &input, const string &input_tensor, const double &mJ) { double J=mJ; if (input_tensor == "b") { - if (fabs(mJ) < sim_iota) { - J = sqrt(det(input)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(input)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside pstretch."); + } + } + bool success_eig_sym = eig_sym(lambda, n_pvectors, input); + if (!success_eig_sym) { + throw simcoon::exception_eig_sym("Error in eig_sym function inside pstretch."); } - eig_sym(lambda, n_pvector, input); lambda.transform( [](double val) { return (sqrt(val)); } ); } else if (input_tensor == "v" || input_tensor == "V") { - if (fabs(mJ) < sim_iota) { - J = det(input); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(input)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside pstretch."); + } + } + bool success_eig_sym = eig_sym(lambda, n_pvectors, input); + if (!success_eig_sym) { + throw simcoon::exception_eig_sym("Error in eig_sym function inside pstretch."); } - eig_sym(lambda, n_pvector, input); } else { throw std::invalid_argument("Invalid input string to describe the input vector: it should be *b* for left Cauchy-Green tensor or *v* or *V* for Eulerian stretch tensor"); } } -void pstretch(vec &lambda, mat &n_pvector, std::vector &N_projectors, const mat &input, const string &input_tensor, const double &mJ) { +void pstretch(vec &lambda, mat &n_pvectors, std::vector &N_projectors, const mat &input, const string &input_tensor, const double &mJ) { - pstretch(lambda, n_pvector, input, input_tensor, mJ); - N_projectors[0] = n_pvector.col(0)*(n_pvector.col(0)).t(); - N_projectors[1] = n_pvector.col(1)*(n_pvector.col(1)).t(); - N_projectors[2] = n_pvector.col(2)*(n_pvector.col(2)).t(); + pstretch(lambda, n_pvectors, input, input_tensor, mJ); + N_projectors[0] = n_pvectors.col(0)*(n_pvectors.col(0)).t(); + N_projectors[1] = n_pvectors.col(1)*(n_pvectors.col(1)).t(); + N_projectors[2] = n_pvectors.col(2)*(n_pvectors.col(2)).t(); } void isochoric_pstretch(vec &lambda_bar, mat &n_pvectors, const mat &input, const string &input_tensor, const double &mJ) { @@ -125,17 +175,33 @@ void isochoric_pstretch(vec &lambda_bar, mat &n_pvectors, const mat &input, cons vec lambda = zeros(3); if (input_tensor == "b") { - if (fabs(mJ) < sim_iota) { - J = sqrt(det(input)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(input)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside isochoric_pstretch."); + } + } + bool success_eig_sym = eig_sym(lambda, n_pvectors, input); + if (!success_eig_sym) { + throw simcoon::exception_eig_sym("Error in eig_sym function inside isochoric_pstretch."); } - eig_sym(lambda, n_pvectors, input); lambda.transform( [](double val) { return (sqrt(val)); } ); } else if (input_tensor == "v" || input_tensor == "V") { - if (fabs(mJ) < sim_iota) { - J = det(input); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(input)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside isochoric_pstretch."); + } } - eig_sym(lambda, n_pvectors, input); + bool success_eig_sym = eig_sym(lambda, n_pvectors, input); + if (!success_eig_sym) { + throw simcoon::exception_eig_sym("Error in eig_sym function inside isochoric_pstretch."); + } } else { throw std::invalid_argument("Invalid input string to describe the input vector: it should be *b* for left Cauchy-Green tensor or *v* or *V* for Eulerian stretch tensor"); @@ -203,8 +269,13 @@ mat tau_iso_hyper_pstretch(const vec &dWdlambda_bar, const vec &lambda_bar, cons mat sigma_iso_hyper_pstretch(const vec &dWdlambda_bar, const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = sqrt(det(b)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(b)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside sigma_iso_hyper_pstretch."); + } } vec lambda_bar = zeros(3); vec n_pvectors = zeros(3); @@ -264,37 +335,73 @@ vec delta_coefs(const vec &a_coefs, const vec &b_coefs, const mat &b) { mat tau_iso_hyper_invariants(const double &dWdI_1_bar, const double &dWdI_2_bar, const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = sqrt(det(b)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(b)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside tau_iso_hyper_invariants."); + } } mat b_bar = pow(J,-2./3.)*b; - mat b_bar2 = powmat(b_bar,2); + + mat b_bar2; + try { + b_bar2 = powmat(b_bar,2); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_powmat("Error in powmat function inside tau_iso_hyper_invariants."); + } + return 2.*dWdI_1_bar*dev(b_bar) + 2.*dWdI_2_bar*(trace(b_bar)*dev(b_bar) - dev(b_bar2)); } mat tau_vol_hyper(const double &dUdJ, const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = sqrt(det(b)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(b)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside tau_vol_hyper."); + } } mat Id = eye(3,3); + return J*dUdJ*eye(3,3); } mat sigma_iso_hyper_invariants(const double &dWdI_1_bar, const double &dWdI_2_bar, const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = sqrt(det(b)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(b)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside sigma_iso_hyper_invariants."); + } } mat b_bar = pow(J,-2./3.)*b; - mat b_bar2 = powmat(b_bar,2); + mat b_bar2; + try { + b_bar2 = powmat(b_bar,2); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_powmat("Error in powmat function inside sigma_iso_hyper_invariants."); + } + return (1./J)*(2.*dWdI_1_bar*dev(b_bar) + 2.*dWdI_2_bar*(trace(b_bar)*dev(b_bar) - dev(b_bar2))); } mat sigma_vol_hyper(const double &dUdJ, const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = sqrt(det(b)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(b)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside sigma_vol_hyper."); + } } mat Id = eye(3,3); return dUdJ*eye(3,3); @@ -302,7 +409,7 @@ mat sigma_vol_hyper(const double &dUdJ, const mat &b, const double &mJ) { /*mat L_iso_hyper_invariants(const double &dWdI_1_bar, const double &dWdI_2_bar, const double &dW2dI_11_bar, const double &dW2dI_12_bar, const double &dW2dI_22_bar, const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { + if (fabs(mJ) < simcoon::iota) { J = sqrt(det(b)); } mat b_bar = pow(J,-2./3.)*b; @@ -328,7 +435,7 @@ mat sigma_vol_hyper(const double &dUdJ, const mat &b, const double &mJ) { /*mat L_iso_hyper_invariants(const vec &delta_coefs, const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { + if (fabs(mJ) < simcoon::iota) { J = sqrt(det(b)); } mat b_bar = pow(J,-2./3.)*b; @@ -382,10 +489,22 @@ mat L_iso_hyper_pstretch(const vec &dWdlambda_bar, const mat &dW2dlambda_bar2, c mat L_iso_hyper_pstretch(const vec &dWdlambda_bar, const mat &dW2dlambda_bar2, const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = sqrt(det(b)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(b)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside L_iso_hyper_pstretch."); + } + } + + vec lambda; + try { + lambda = eig_sym(b); + } catch (const std::runtime_error &e) { + cerr << "Error in eig_sym: " << e.what() << endl; + throw simcoon::exception_eig_sym("Failed to compute eigenvalues in L_iso_hyper_pstretch."); } - vec lambda = eig_sym(b); vec lambda_bar = zeros(3); vec n_pvectors = zeros(3); @@ -420,11 +539,24 @@ mat L_iso_hyper_pstretch(const vec &dWdlambda_bar, const mat &dW2dlambda_bar2, c mat L_iso_hyper_invariants(const double &dWdI_1_bar, const double &dWdI_2_bar, const double &dW2dI_11_bar, const double &dW2dI_12_bar, const double &dW2dI_22_bar, const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = sqrt(det(b)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(b)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside L_iso_hyper_invariants."); + } } mat b_bar = pow(J,-2./3.)*b; - mat b_bar2 = powmat(b_bar,2); + + mat b_bar2; + try { + b_bar2 = powmat(b_bar,2); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_powmat("Error in powmat function inside L_iso_hyper_invariants."); + } + mat dev_b_bar = dev(b_bar); mat dev_b_bar2 = dev(b_bar2); vec I_bar = isochoric_invariants(b,J); @@ -445,8 +577,13 @@ mat L_iso_hyper_invariants(const double &dWdI_1_bar, const double &dWdI_2_bar, c mat L_vol_hyper(const double &dUdJ, const double &dU2dJ2, const mat &b, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = sqrt(det(b)); + if (fabs(mJ) < simcoon::iota) { + try { + J = sqrt(det(b)); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside L_vol_hyper."); + } } return (dUdJ+dU2dJ2*J)*3.*Ivol() - 2.*dUdJ*Ireal(); } diff --git a/src/Continuum_mechanics/Functions/kinematics.cpp b/src/Continuum_mechanics/Functions/kinematics.cpp index 3a638ad95..64264433b 100755 --- a/src/Continuum_mechanics/Functions/kinematics.cpp +++ b/src/Continuum_mechanics/Functions/kinematics.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include using namespace std; @@ -49,7 +50,13 @@ mat ER_to_F(const mat &E, const mat &R) { vec N = N_alpha.col(i); U = U + (lambda_alpha(i)*(N*N.t())); }*/ - return (R*sqrtmat_sympd(C)); + + try { + return (R*sqrtmat_sympd(C)); + } catch (const std::runtime_error &e) { + cerr << "Error in sqrtmat_sympd : " << e.what() << endl; + throw simcoon::exception_sqrtmat_sympd("Error in sqrtmat_sympd function inside ER_to_F."); + } } mat eR_to_F(const mat &e, const mat &R) { @@ -60,7 +67,13 @@ mat eR_to_F(const mat &e, const mat &R) { assert(R.n_rows == 3); //From e we compute V : e = 1/2 (C-I) --> C = U^2 = 2E+I - mat V = expmat_sym(e); + mat V; + try { + V = expmat_sym(e); + } catch (const std::runtime_error &e) { + cerr << "Error in expmat_sym : " << e.what() << endl; + throw simcoon::exception_expmat_sym("Error in expmat_sym function inside eR_to_F."); + } return (V*R); } @@ -72,7 +85,13 @@ mat G_UdX(const mat &F) { //This function computes the gradient of displacement (Eulerian) from the deformation gradient tensor mat G_Udx(const mat &F) { - return eye(3,3) - inv(F); + + try { + return eye(3,3) - inv(F); + } catch (const std::runtime_error &e) { + cerr << "Error in inv : " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside G_Udx."); + } } //This function computes the Right Cauchy-Green C @@ -88,15 +107,39 @@ mat L_Cauchy_Green(const mat &F) { //Provides the RU decomposition of the transformation gradient F void RU_decomposition(mat &R, mat &U, const mat &F) { mat U2 = F.t()*F; - U = sqrtmat_sympd(U2); - R = F*inv(U); + + try { + U = sqrtmat_sympd(U2); + } catch (const std::runtime_error &e) { + cerr << "Error in sqrtmat_sympd : " << e.what() << endl; + throw simcoon::exception_sqrtmat_sympd("Error in sqrtmat_sympd function inside RU_decomposition."); + } + + try { + R = F*inv(U); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside RU_decomposition."); + } } //Provides the VR decomposition of the transformation gradient F void VR_decomposition(mat &V, mat &R, const mat &F) { mat V2 = F*F.t(); - V = sqrtmat_sympd(V2); - R = inv(V)*F; + + try { + V = sqrtmat_sympd(V2); + } catch (const std::runtime_error &e) { + cerr << "Error in sqrtmat_sympd : " << e.what() << endl; + throw simcoon::exception_sqrtmat_sympd("Error in sqrtmat_sympd function inside VR_decomposition."); + } + + try { + R = inv(V)*F; + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside RU_decomposition."); + } } //This function computes the common Right (or Left) Cauchy-Green invariants @@ -105,13 +148,18 @@ vec Inv_X(const mat &X) { vec I = zeros(3); I(0) = trace(X); I(1) = 0.5*(pow(trace(X),2.) - trace(X*X)); - I(2) = det(X); + I(2) = det(X); return I; } //This function computes the Cauchy deformation tensor b mat Cauchy(const mat &F) { - return inv(L_Cauchy_Green(F)); + try { + return inv(L_Cauchy_Green(F)); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Cauchy."); + } } //This function computes the Green-Lagrange finite strain tensor E @@ -126,22 +174,37 @@ mat Euler_Almansi(const mat &F) { //This function computes the Euler-Almansi finite strain tensor h mat Log_strain(const mat &F) { - return 0.5*logmat_sympd(L_Cauchy_Green(F)); + try { + return 0.5*logmat_sympd(L_Cauchy_Green(F)); + } catch (const std::runtime_error &e) { + cerr << "Error in logmat_sympd: " << e.what() << endl; + throw simcoon::exception_logmat_sympd("Error in logmat_sympd function inside Log_strain."); + } } //This function computes the velocity difference mat finite_L(const mat &F0, const mat &F1, const double &DTime) { //Definition of L = dot(F)*F^-1 - return (1./DTime)*(F1-F0)*inv(F1); + try { + return (1./DTime)*(F1-F0)*inv(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside finite_L."); + } } //This function computes the deformation rate D mat finite_D(const mat &F0, const mat &F1, const double &DTime) { //Definition of L = dot(F)*F^-1 - mat L = (1./DTime)*(F1-F0)*inv(F1); - + mat L; + try { + L = (1./DTime)*(F1-F0)*inv(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside finite_D."); + } //Definition of the deformation rate D return 0.5*(L+L.t()); @@ -151,7 +214,13 @@ mat finite_D(const mat &F0, const mat &F1, const double &DTime) { mat finite_W(const mat &F0, const mat &F1, const double &DTime) { //Definition of L = dot(F)*F^-1 - mat L = (1./DTime)*(F1-F0)*inv(F1); + mat L; + try { + L = (1./DTime)*(F1-F0)*inv(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside finite_W."); + } //Definition of the rotation matrix Q return 0.5*(L-L.t()); @@ -175,7 +244,13 @@ mat finite_Omega(const mat &F0, const mat &F1, const double &DTime) { //This function computes the increment of finite rotation mat finite_DQ(const mat &Omega0, const mat &Omega1, const double &DTime) { - return (eye(3,3)+0.5*DTime*Omega0)*(inv(eye(3,3)-0.5*DTime*Omega1)); + try { + return (eye(3,3)+0.5*DTime*Omega0)*(inv(eye(3,3)-0.5*DTime*Omega1)); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside finite_DQ."); + } + } diff --git a/src/Continuum_mechanics/Functions/natural_basis.cpp b/src/Continuum_mechanics/Functions/natural_basis.cpp index 314d06b93..7b303b43a 100755 --- a/src/Continuum_mechanics/Functions/natural_basis.cpp +++ b/src/Continuum_mechanics/Functions/natural_basis.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include using namespace std; @@ -70,7 +71,12 @@ natural_basis::natural_basis(const std::vector &mg_i) } } - g0ij = inv(g_ij); + try { + g0ij = inv(g_ij); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside natural_basis constructor."); + } for (unsigned int i=0; i < 3; i++) { for (unsigned int j=0; j<3; j++) { @@ -130,7 +136,12 @@ void natural_basis::update(const std::vector &mg_i) } } - g0ij = inv(g_ij); + try { + g0ij = inv(g_ij); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside natural_basis update."); + } for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { @@ -158,7 +169,12 @@ void natural_basis::from_F(const arma::mat &F) } } - g0ij = inv(g_ij); + try { + g0ij = inv(g_ij); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside natural_basis from_F."); + } for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { diff --git a/src/Continuum_mechanics/Functions/objective_rates.cpp b/src/Continuum_mechanics/Functions/objective_rates.cpp index 65bcb5bcf..a57082636 100755 --- a/src/Continuum_mechanics/Functions/objective_rates.cpp +++ b/src/Continuum_mechanics/Functions/objective_rates.cpp @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -39,14 +40,28 @@ namespace simcoon{ void Jaumann(mat &DR, mat &D, mat &W, const double &DTime, const mat &F0, const mat &F1) { mat I = eye(3,3); - mat L = (1./DTime)*(F1-F0)*inv(F1); + mat L; + if(DTime > simcoon::iota) { + try { + L = (1./DTime)*(F1-F0)*inv(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Jaumann (L)."); + } + } + //decomposition of L D = 0.5*(L+L.t()); W = 0.5*(L-L.t()); //Jaumann - DR = (inv(I-0.5*DTime*W))*(I+0.5*DTime*W); + try { + DR = (inv(I-0.5*DTime*W))*(I+0.5*DTime*W); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Jaumann (DR)."); + } } void Green_Naghdi(mat &DR, mat &D, mat &Omega, const double &DTime, const mat &F0, const mat &F1) { @@ -59,14 +74,28 @@ void Green_Naghdi(mat &DR, mat &D, mat &Omega, const double &DTime, const mat &F RU_decomposition(R0,U0,F0); RU_decomposition(R1,U1,F1); - mat L = (1./DTime)*(F1-F0)*inv(F1); + mat L; + if(DTime > simcoon::iota) { + try { + L = (1./DTime)*(F1-F0)*inv(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Green_Naghdi (L)."); + } + } //decomposition of L D = 0.5*(L+L.t()); mat W = 0.5*(L-L.t()); Omega = (1./DTime)*(R1-R0)*R1.t(); - DR = (inv(I-0.5*DTime*Omega))*(I+0.5*DTime*Omega); + + 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 Green_Naghdi (DR)."); + } //alternative ... to test // DR = (F1-F0)*inv(U1)-R0*(U1-U0)*inv(U1); } @@ -81,14 +110,27 @@ void logarithmic_R(mat &DR, mat &N_1, mat &N_2, mat &D, mat &Omega, const double RU_decomposition(R0,U0,F0); RU_decomposition(R1,U1,F1); - mat L = (1./DTime)*(F1-F0)*inv(F1); + mat L; + if(DTime > simcoon::iota) { + try { + L = (1./DTime)*(F1-F0)*inv(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside logarithmic_R (L)."); + } + } //decomposition of L D = 0.5*(L+L.t()); mat W = 0.5*(L-L.t()); Omega = (1./DTime)*(R1-R0)*R1.t(); - DR = (inv(I-0.5*DTime*Omega))*(I+0.5*DTime*Omega); + 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 logarithmic_R (DR)."); + } //alternative ... to test // DR = (F1-F0)*inv(U1)-R0*(U1-U0)*inv(U1); @@ -97,7 +139,11 @@ void logarithmic_R(mat &DR, mat &N_1, mat &N_2, mat &D, mat &Omega, const double vec bi = zeros(3); mat Bi; - eig_sym(bi, Bi, B); + bool success_eig_sym = eig_sym(bi, Bi, B); + if (!success_eig_sym) { + throw simcoon::exception_eig_sym("Error in eig_sym function inside logarithmic_R."); + } + std::vector Bi_proj(3); Bi_proj[0] = Bi.col(0)*(Bi.col(0)).t(); Bi_proj[1] = Bi.col(1)*(Bi.col(1)).t(); @@ -106,7 +152,7 @@ void logarithmic_R(mat &DR, mat &N_1, mat &N_2, mat &D, mat &Omega, const double N_1 = zeros(3,3); for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { - if ((i!=j)&&(fabs(bi(i)-bi(j))>sim_iota)) { + if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) { N_1+=((1.+(bi(i)/bi(j)))/(1.-(bi(i)/bi(j)))+2./log(bi(i)/bi(j)))*Bi_proj[i]*D*Bi_proj[j]; } } @@ -115,7 +161,7 @@ void logarithmic_R(mat &DR, mat &N_1, mat &N_2, mat &D, mat &Omega, const double N_2 = zeros(3,3); for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { - if ((i!=j)&&(fabs(bi(i)-bi(j))>sim_iota)) { + if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) { N_2+=((1.-(pow(bi(i)/bi(j),0.5)))/(1.+(pow(bi(i)/bi(j),0.5))))*Bi_proj[i]*D*Bi_proj[j]; } } @@ -132,20 +178,28 @@ void logarithmic_F(mat &DF, mat &N_1, mat &N_2, mat &D, mat &L, const double &DT RU_decomposition(R0,U0,F0); RU_decomposition(R1,U1,F1); - L = (1./DTime)*(F1-F0)*inv(F1); + if(DTime > simcoon::iota) { + try { + L = (1./DTime)*(F1-F0)*inv(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside logarithmic_F (L)."); + } + } //decomposition of L D = 0.5*(L+L.t()); mat W = 0.5*(L-L.t()); - - //alternative ... to test //Logarithmic mat B = L_Cauchy_Green(F1); vec bi = zeros(3); mat Bi; - eig_sym(bi, Bi, B); + bool success_eig_sym = eig_sym(bi, Bi, B); + if (!success_eig_sym) { + throw simcoon::exception_eig_sym("Error in eig_sym function inside logarithmic_R."); + } std::vector Bi_proj(3); Bi_proj[0] = Bi.col(0)*(Bi.col(0)).t(); Bi_proj[1] = Bi.col(1)*(Bi.col(1)).t(); @@ -154,7 +208,7 @@ void logarithmic_F(mat &DF, mat &N_1, mat &N_2, mat &D, mat &L, const double &DT N_1 = zeros(3,3); for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { - if ((i!=j)&&(fabs(bi(i)-bi(j))>sim_iota)) { + if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) { N_1+=((1.+(bi(i)/bi(j)))/(1.-(bi(i)/bi(j)))+2./log(bi(i)/bi(j)))*Bi_proj[i]*D*Bi_proj[j]; } } @@ -163,24 +217,41 @@ void logarithmic_F(mat &DF, mat &N_1, mat &N_2, mat &D, mat &L, const double &DT N_2 = zeros(3,3); for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { - if ((i!=j)&&(fabs(bi(i)-bi(j))>sim_iota)) { + if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) { N_2+=((1.-(pow(bi(i)/bi(j),0.5)))/(1.+(pow(bi(i)/bi(j),0.5))))*Bi_proj[i]*D*Bi_proj[j]; } } } - DF = (inv(I-0.5*DTime*L))*(I+0.5*DTime*L); - + try { + DF = (inv(I-0.5*DTime*L))*(I+0.5*DTime*L); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside logarithmic_F (DF)."); + } } void Truesdell(mat &DF, mat &D, mat &L, const double &DTime, const mat &F0, const mat &F1) { mat I = eye(3,3); - L = (1./DTime)*(F1-F0)*inv(F1); + if(DTime > simcoon::iota) { + try { + L = (1./DTime)*(F1-F0)*inv(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Truesdell (L)."); + } + } + //Note that The "spin" is actually L (spin for rigid frames of reference, "flot" for Truesdell) D = 0.5*(L+L.t()); //Truesdell - DF = (inv(I-0.5*DTime*L))*(I+0.5*DTime*L); + try { + DF = (inv(I-0.5*DTime*L))*(I+0.5*DTime*L); + } 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 get_BBBB(const mat &F1) { @@ -188,13 +259,16 @@ mat get_BBBB(const mat &F1) { vec bi = zeros(3); mat Bi; - eig_sym(bi, Bi, B); + bool success_eig_sym = eig_sym(bi, Bi, B); + if (!success_eig_sym) { + throw simcoon::exception_eig_sym("Error in eig_sym function inside logarithmic_R."); + } mat BBBB = zeros(6,6); double f_z = 0.; for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { - if ((i!=j)&&(fabs(bi(i)-bi(j))>sim_iota)) { + if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) { f_z = (1.+(bi(i)/bi(j)))/(1.-(bi(i)/bi(j)))+2./log(bi(i)/bi(j)); BBBB = BBBB + f_z*B_klmn(Bi.col(i),Bi.col(j)); } @@ -208,13 +282,16 @@ mat get_BBBB_GN(const mat &F1) { vec bi = zeros(3); mat Bi; - eig_sym(bi, Bi, B); + bool success_eig_sym = eig_sym(bi, Bi, B); + if (!success_eig_sym) { + throw simcoon::exception_eig_sym("Error in eig_sym function inside logarithmic_R."); + } mat BBBB = zeros(6,6); double f_z = 0.; for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { - if ((i!=j)&&(fabs(bi(i)-bi(j))>sim_iota)) { + if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) { f_z = (sqrt(bi(j)) - sqrt(bi(i)))/(sqrt(bi(j)) + sqrt(bi(i))); BBBB = BBBB + f_z*B_klmn(Bi.col(i),Bi.col(j)); } @@ -226,10 +303,15 @@ mat get_BBBB_GN(const mat &F1) { void logarithmic(mat &DR, mat &D, mat &Omega, const double &DTime, const mat &F0, const mat &F1) { mat I = eye(3,3); mat L = zeros(3,3); - - if(DTime > sim_iota) { - L = (1./DTime)*(F1-F0)*inv(F1); - } + + if(DTime > simcoon::iota) { + try { + L = (1./DTime)*(F1-F0)*inv(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside logarithmic (L)."); + } + } //decomposition of L D = 0.5*(L+L.t()); @@ -240,7 +322,10 @@ void logarithmic(mat &DR, mat &D, mat &Omega, const double &DTime, const mat &F0 vec bi = zeros(3); mat Bi; - eig_sym(bi, Bi, B); + bool success_eig_sym = eig_sym(bi, Bi, B); + if (!success_eig_sym) { + throw simcoon::exception_eig_sym("Error in eig_sym function inside logarithmic_R."); + } std::vector Bi_proj(3); Bi_proj[0] = Bi.col(0)*(Bi.col(0)).t(); Bi_proj[1] = Bi.col(1)*(Bi.col(1)).t(); @@ -249,25 +334,43 @@ void logarithmic(mat &DR, mat &D, mat &Omega, const double &DTime, const mat &F0 mat N = zeros(3,3); for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { - if ((i!=j)&&(fabs(bi(i)-bi(j))>sim_iota)) { + if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) { N+=((1.+(bi(i)/bi(j)))/(1.-(bi(i)/bi(j)))+2./log(bi(i)/bi(j)))*Bi_proj[i]*D*Bi_proj[j]; } } } Omega = W + N; - DR = (inv(I-0.5*DTime*Omega))*(I+0.5*DTime*Omega); + + 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 logarithmic (DR)."); + } } mat Delta_log_strain(const mat &D, const mat &Omega, const double &DTime) { mat I = eye(3,3); - mat DR = (inv(I-0.5*DTime*Omega))*(I+0.5*DTime*Omega); + 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 logarithmic (DR)."); + } return 0.5*(D+(DR*D*DR.t()))*DTime; } //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 mat DtauDe_2_DSDE(const mat &Lt, const mat &B, const mat &F, const mat &tau){ - mat invF = inv(F); + mat invF; + try { + invF = inv(F); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside DtauDe_2_DSDE."); + } Tensor2 invF_ = mat_FTensor2(invF); Tensor2 delta_ = mat_FTensor2(eye(3,3)); Tensor2 tau_ = mat_FTensor2(tau); @@ -296,7 +399,14 @@ mat DtauDe_2_DSDE(const mat &Lt, const mat &B, const mat &F, const mat &tau){ mat Dtau_LieDD_2_DSDE(const mat &Lt, const mat &F){ - mat invF = inv(F); + mat invF; + try { + invF = inv(F); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Dtau_LieDD_2_DSDE."); + } + Tensor2 invF_ = mat_FTensor2(invF); Tensor4 Dtau_LieDD_ = mat_FTensor4(Lt); Tensor4 DSDE_ = mat_FTensor4(zeros(6,6)); @@ -317,7 +427,13 @@ mat Dtau_LieDD_2_DSDE(const mat &Lt, const mat &F){ mat DtauDe_JaumannDD_2_DSDE(const mat &Lt, const mat &F, const mat &tau){ - mat invF = inv(F); + mat invF; + try { + invF = inv(F); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside DtauDe_JaumannDD_2_DSDE."); + } Tensor2 invF_ = mat_FTensor2(invF); Tensor2 delta_ = mat_FTensor2(eye(3,3)); Tensor2 tau_ = mat_FTensor2(tau); @@ -345,38 +461,64 @@ mat DtauDe_JaumannDD_2_DSDE(const mat &Lt, const mat &F, const mat &tau){ //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 mat DsigmaDe_2_DSDE(const mat &Lt, const mat &B, const mat &F, const mat &sigma){ - double J = det(F); + 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_2_DSDE."); + } return J*DtauDe_2_DSDE(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 mat DsigmaDe_2_DSDE(const mat &Lt, const mat &F, const mat &sigma){ - double J = det(F); + 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_2_DSDE."); + } mat B = get_BBBB(F); return J*DtauDe_2_DSDE(Lt, B, F, Cauchy2Kirchoff(sigma, F, J)); } mat Dsigma_LieDD_2_DSDE(const mat &Lt, const mat &F){ - double J = det(F); + 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); } mat DsigmaDe_JaumannDD_2_DSDE(const mat &Lt, const mat &F, const mat &sigma){ - double J = det(F); + 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)); } mat DtauDe_2_DsigmaDe(const mat &Lt, const double &J) { + assert(J > simcoon::iota); return (1./J)*Lt; } mat DsigmaDe_2_DtauDe(const mat &Lt, const double &J) { + assert(J > simcoon::iota); return Lt*J; } @@ -408,7 +550,13 @@ mat DSDE_2_DtauDe(const mat &DSDE, const mat &B, const mat &F, const mat &tau) { mat DSDE_2_DsigmaDe(const mat &DSDE, const mat &B, const mat &F, const mat &sigma) { - double J = det(F); + 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_DsigmaDe."); + } return (1./J)*DSDE_2_DtauDe(DSDE, B, F, Cauchy2Kirchoff(sigma, F, J)); } @@ -435,7 +583,13 @@ mat DSDE_2_Dtau_LieDD(const mat &DSDE, const mat &F) { mat DSDE_2_DsigmaDe_LieDD(const mat &DSDE, const mat &F) { - double J = det(F); + 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_DsigmaDe_LieDD."); + } return (1./J)*DSDE_2_Dtau_LieDD(DSDE, F); } @@ -464,7 +618,13 @@ 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 = det(F); + 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)); } diff --git a/src/Continuum_mechanics/Functions/recovery_props.cpp b/src/Continuum_mechanics/Functions/recovery_props.cpp index 4c5f6a0ea..c77041af4 100755 --- a/src/Continuum_mechanics/Functions/recovery_props.cpp +++ b/src/Continuum_mechanics/Functions/recovery_props.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -34,9 +35,9 @@ namespace simcoon{ void check_symetries(const mat &L, std::string &umat_type, int &axis, vec &props, int &maj_sym, const double &tol) { - double max_tol_sim_limit = sim_limit; - if(tol > sim_limit) { - max_tol_sim_limit = tol; + double max_tol_sim = simcoon::limit; + if(tol > simcoon::limit) { + max_tol_sim = tol; } axis = 0; //Indicate no preferential axis @@ -92,14 +93,14 @@ void check_symetries(const mat &L, std::string &umat_type, int &axis, vec &props for (unsigned int i=0; i #include #include +#include #include #include #include @@ -33,15 +34,25 @@ namespace simcoon{ mat Cauchy2PKI(const mat &sigma, const mat &F, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = det(F); + if (fabs(mJ) < simcoon::iota) { + 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 Cauchy2PKI."); + } } //If J is still less than a small value, we assume that sigma=tau=PK1=PKII = 0 - if (fabs(J) < sim_iota) { + if (fabs(J) < simcoon::iota) { return zeros(3,3); } else { - return J*sigma*inv(F.t()); + try { + return J*sigma*inv(F.t()); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Cauchy2PKI."); + } } } @@ -51,7 +62,7 @@ mat Cauchy2Biot(const mat &sigma, const mat &F, const mat &mR, const double &mJ) mat PKI = Cauchy2PKI(sigma, F, mJ); mat R = mR; - if (norm(R,1) < sim_iota) { + if (norm(R,1) < simcoon::iota) { mat U = zeros(3,3); RU_decomposition(R,U,F); } @@ -62,15 +73,25 @@ mat Cauchy2Biot(const mat &sigma, const mat &F, const mat &mR, const double &mJ) mat Cauchy2PKII(const mat &sigma, const mat &F, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = det(F); + if (fabs(mJ) < simcoon::iota) { + 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 Cauchy2PKII."); + } } //If J is still less than a small value, we assume that sigma=tau=PK1=PKII = 0 - if (fabs(J) < sim_iota) { + if (fabs(J) < simcoon::iota) { return zeros(3,3); } else { - return J*inv(F)*sigma*inv(F.t()); + try { + return J*inv(F)*sigma*inv(F.t()); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Cauchy2PKII."); + } } } @@ -78,11 +99,16 @@ mat Cauchy2PKII(const mat &sigma, const mat &F, const double &mJ) { mat Cauchy2Kirchoff(const mat &sigma, const mat &F, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = det(F); + if (fabs(mJ) < simcoon::iota) { + 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 Cauchy2Kirchoff."); + } } //If J is still less than a small value, we assume that sigma=tau=PK1=PKII = 0 - if (fabs(J) < sim_iota) { + if (fabs(J) < simcoon::iota) { return zeros(3,3); } else { @@ -95,11 +121,16 @@ mat Cauchy2Kirchoff(const mat &sigma, const mat &F, const double &mJ) { vec Cauchy2Kirchoff(const vec &sigma, const mat &F, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = det(F); + if (fabs(mJ) < simcoon::iota) { + 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 Cauchy2Kirchoff."); + } } //If J is still less than a small value, we assume that sigma=tau=PK1=PKII = 0 - if (fabs(J) < sim_iota) { + if (fabs(J) < simcoon::iota) { return zeros(6); } else { @@ -112,11 +143,16 @@ vec Cauchy2Kirchoff(const vec &sigma, const mat &F, const double &mJ) { mat Kirchoff2Cauchy(const mat &tau, const mat &F, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = det(F); + if (fabs(mJ) < simcoon::iota) { + 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 Kirchoff2Cauchy."); + } } //If J is still less than a small value, we assume that sigma=tau=PK1=PKII = 0 - if (fabs(J) < sim_iota) { + if (fabs(J) < simcoon::iota) { return zeros(3,3); } else { @@ -128,11 +164,16 @@ mat Kirchoff2Cauchy(const mat &tau, const mat &F, const double &mJ) { vec Kirchoff2Cauchy(const vec& tau, const mat& F, const double& mJ) { double J = mJ; - if (fabs(mJ) < sim_iota) { - J = det(F); + if (fabs(mJ) < simcoon::iota) { + 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 Kirchoff2Cauchy."); + } } //If J is still less than a small value, we assume that sigma=tau=PK1=PKII = 0 - if (fabs(J) < sim_iota) { + if (fabs(J) < simcoon::iota) { return zeros(6); } else { @@ -145,20 +186,36 @@ vec Kirchoff2Cauchy(const vec& tau, const mat& F, const double& mJ) { mat Kirchoff2PKI(const mat &tau, const mat &F, const double &mJ) { UNUSED(mJ); - return tau*inv(F.t()); + try { + return tau*inv(F.t()); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Kirchoff2PKI."); + } } //This function returns the second Piola-Kirchoff stress tensor from the Kirchoff stress tensor, the transformation gradient F and its determinant (optional, if not indicated, it will be computed) mat Kirchoff2PKII(const mat &tau, const mat &F, const double &mJ) { UNUSED(mJ); - return inv(F)*tau*inv(F.t()); + try { + return inv(F)*tau*inv(F.t()); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Kirchoff2PKII."); + } + } vec Kirchoff2PKII(const vec &tau, const mat &F, const double &mJ) { UNUSED(mJ); - return t2v_stress(inv(F)*v2t_stress(tau)*inv(F.t())); + try { + return t2v_stress(inv(F)*v2t_stress(tau)*inv(F.t())); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Kirchoff2PKII."); + } } @@ -180,11 +237,16 @@ mat PKII2Kirchoff(const mat &S, const mat &F, const double &mJ) { mat PKI2Cauchy(const mat &Sigma, const mat &F, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = det(F); + if (fabs(mJ) < simcoon::iota) { + 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 PKI2Cauchy."); + } } //If J is still less than a small value, we assume that sigma=tau=PK1=PKII = 0 - if (fabs(J) < sim_iota) { + if (fabs(J) < simcoon::iota) { return zeros(3,3); } else { @@ -196,11 +258,16 @@ mat PKI2Cauchy(const mat &Sigma, const mat &F, const double &mJ) { mat PKII2Cauchy(const mat &S, const mat &F, const double &mJ) { double J=mJ; - if (fabs(mJ) < sim_iota) { - J = det(F); + if (fabs(mJ) < simcoon::iota) { + 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 PKII2Cauchy."); + } } //If J is still less than a small value, we assume that sigma=tau=PK1=PKII = 0 - if (fabs(J) < sim_iota) { + if (fabs(J) < simcoon::iota) { return zeros(3,3); } else { diff --git a/src/Continuum_mechanics/Homogenization/ellipsoid_multi.cpp b/src/Continuum_mechanics/Homogenization/ellipsoid_multi.cpp index 56581ee89..d8d2fdc3f 100755 --- a/src/Continuum_mechanics/Homogenization/ellipsoid_multi.cpp +++ b/src/Continuum_mechanics/Homogenization/ellipsoid_multi.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -133,8 +134,12 @@ void ellipsoid_multi::fillT(const mat& Lt_m, const mat& Lt, const ellipsoid &ell S_loc = Eshelby(Lt_m_local_geom, ell.a1, ell.a2, ell.a3, x, wx, y, wy, mp, np); mat Lt_local_geom = rotate_g2l_L(Lt, ell.psi_geom, ell.theta_geom, ell.phi_geom); - T_loc = inv(eye(6,6) + S_loc*inv(Lt_m_local_geom)*(Lt_local_geom - Lt_m_local_geom)); - + try { + T_loc = inv(eye(6,6) + S_loc*inv(Lt_m_local_geom)*(Lt_local_geom - Lt_m_local_geom)); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside ellipsoid_multi::fillT."); + } T = rotate_l2g_A(T_loc, ell.psi_geom, ell.theta_geom, ell.phi_geom); } @@ -147,7 +152,12 @@ void ellipsoid_multi::fillT_iso(const mat& Lt_m, const mat& Lt, const ellipsoid S_loc = Eshelby(Lt_m_iso, ell.a1, ell.a2, ell.a3, x, wx, y, wy, mp, np); mat Lt_local_geom = rotate_g2l_L(Lt, ell.psi_geom, ell.theta_geom, ell.phi_geom); - T_loc = inv(eye(6,6) + S_loc*inv(Lt_m_iso)*(Lt_local_geom - Lt_m_iso)); + try { + T_loc = inv(eye(6,6) + S_loc*inv(Lt_m_iso)*(Lt_local_geom - Lt_m_iso)); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside ellipsoid_multi::fillT_is."); + } T = rotate_l2g_A(T_loc, ell.psi_geom, ell.theta_geom, ell.phi_geom); } @@ -162,10 +172,15 @@ void ellipsoid_multi::fillT_mec_in(const mat& L_m, const mat& L, const ellipsoid S_loc = Eshelby(L_m_local_geom, ell.a1, ell.a2, ell.a3, x, wx, y, wy, mp, np); mat L_local_geom = rotate_g2l_L(L, ell.psi_geom, ell.theta_geom, ell.phi_geom); - T_loc = inv(eye(6,6) + S_loc*inv(L_m_local_geom)*(L_local_geom - L_m_local_geom)); + try { + T_loc = inv(eye(6,6) + S_loc*inv(L_m_local_geom)*(L_local_geom - L_m_local_geom)); + T_in_loc = (eye(6,6)-T_loc)*inv(L_m_local_geom - L_local_geom); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside ellipsoid_multi::fillT_mec_in."); + } + T = rotate_l2g_A(T_loc, ell.psi_geom, ell.theta_geom, ell.phi_geom); - - T_in_loc = (eye(6,6)-T_loc)*inv(L_m_local_geom - L_local_geom); T_in = rotate_l2g_M(T_in_loc, ell.psi_geom, ell.theta_geom, ell.phi_geom); } diff --git a/src/Continuum_mechanics/Homogenization/eshelby.cpp b/src/Continuum_mechanics/Homogenization/eshelby.cpp index e01d1eecb..2ceb968b0 100755 --- a/src/Continuum_mechanics/Homogenization/eshelby.cpp +++ b/src/Continuum_mechanics/Homogenization/eshelby.cpp @@ -281,13 +281,13 @@ mat Eshelby(const mat &Lt, const double &a1, const double &a2, const double &a3, for (int i=0; i<3; i++) { for (int j=0; j<6; j++) { - S(i,j) = S(i,j)*(1./(8.*sim_pi)); + S(i,j) = S(i,j)*(1./(8.*simcoon::pi)); } } for (int i=3; i<6; i++) { for (int j=0; j<6; j++) { - S(i,j) = S(i,j)*(1./(4.*sim_pi)); + S(i,j) = S(i,j)*(1./(4.*simcoon::pi)); } } return S; @@ -377,13 +377,13 @@ mat T_II(const mat &Lt, const double &a1, const double &a2, const double &a3, co for (int i=0; i<3; i++) { for (int j=0; j<6; j++) { - T_II(i,j) = T_II(i,j)*(1./(8.*sim_pi)); + T_II(i,j) = T_II(i,j)*(1./(8.*simcoon::pi)); } } for (int i=3; i<6; i++) { for (int j=0; j<6; j++) { - T_II(i,j) = T_II(i,j)*(1./(4.*sim_pi)); + T_II(i,j) = T_II(i,j)*(1./(4.*simcoon::pi)); } } return T_II; @@ -404,7 +404,7 @@ void points(vec &x, vec &wx, vec &y, vec &wy, const int &mp, const int &np) { double x1=0.; - double x2=2.*sim_pi; + double x2=2.*simcoon::pi; int n=0; int m=0; double p1=0.; @@ -433,10 +433,10 @@ void points(vec &x, vec &wx, vec &y, vec &wy, const int &mp, const int &np) //Loop over the desired roots for (int i=1; i<=m; i++) { - z=cos(sim_pi*(i-0.25)/(n+0.5)); + z=cos(simcoon::pi*(i-0.25)/(n+0.5)); z1=0.; - while(fabs(z-z1)>sim_iota) { + while(fabs(z-z1)>simcoon::iota) { p1 = 1.; p2 = 0.; diff --git a/src/Continuum_mechanics/Material/ODF.cpp b/src/Continuum_mechanics/Material/ODF.cpp index 8c777f22c..d66837eab 100755 --- a/src/Continuum_mechanics/Material/ODF.cpp +++ b/src/Continuum_mechanics/Material/ODF.cpp @@ -70,8 +70,8 @@ ODF::ODF(const int &nAngle, const bool &nradian, const double &angle_min, const limits(1) = angle_max; if (radian == false) { - limits(0) *= (sim_pi/180.); - limits(1) *= (sim_pi/180.); + limits(0) *= (simcoon::pi/180.); + limits(1) *= (simcoon::pi/180.); } } diff --git a/src/Continuum_mechanics/Material/ODF2Nphases.cpp b/src/Continuum_mechanics/Material/ODF2Nphases.cpp index 50bc04816..c86a2ca79 100755 --- a/src/Continuum_mechanics/Material/ODF2Nphases.cpp +++ b/src/Continuum_mechanics/Material/ODF2Nphases.cpp @@ -47,12 +47,12 @@ vec get_densities_ODF(const vec &x, const string &path_data, const string &input vec y = zeros(x.n_elem); vec x_rad; if (!radian) { - x_rad = x*(sim_pi/180.); + x_rad = x*(simcoon::pi/180.); if(x_rad.min() < 0.) { cout << "Error : x.min() < 0. Please provide an angle vector with all angles >=0."; return y; } - if(x_rad.max() > sim_pi) { + if(x_rad.max() > simcoon::pi) { cout << "Error : x.max() > pi Please provide an angle vector with all angles <=pi/180"; return y; } @@ -62,7 +62,7 @@ vec get_densities_ODF(const vec &x, const string &path_data, const string &input cout << "Error : x.min() < 0. Please provide an angle vector with all angles >=0."; return y; } - if(x.max() > sim_pi) { + if(x.max() > simcoon::pi) { cout << "Error : x.max() > pi Please provide an angle vector with all angles <=pi"; return y; } @@ -162,7 +162,7 @@ phase_characteristics discretize_ODF(const phase_characteristics &rve_init, ODF fill_angles(alpha, temp, odf_rve, angles_mat); if(alpha - odf_rve.limits(0) < dalpha) - temp.sptr_shape->concentration = dalpha/6. * (odf_rve.density(sim_pi+alpha-dalpha/2) + 4.*odf_rve.density(alpha) + odf_rve.density(alpha+dalpha/2)); + temp.sptr_shape->concentration = dalpha/6. * (odf_rve.density(simcoon::pi+alpha-dalpha/2) + 4.*odf_rve.density(alpha) + odf_rve.density(alpha+dalpha/2)); else temp.sptr_shape->concentration = dalpha/6. * (odf_rve.density(alpha-dalpha/2) + 4.*odf_rve.density(alpha) + odf_rve.density(alpha+dalpha/2)); diff --git a/src/Continuum_mechanics/Material/peak.cpp b/src/Continuum_mechanics/Material/peak.cpp index 16e47cf41..cfd181e59 100755 --- a/src/Continuum_mechanics/Material/peak.cpp +++ b/src/Continuum_mechanics/Material/peak.cpp @@ -117,25 +117,25 @@ double peak::get_density_ODF(const double &theta) break; } case 2: { - return ODF_hard(theta, mean, s_dev, ampl) + ODF_hard(theta - sim_pi, mean, s_dev, ampl) + ODF_hard(theta + sim_pi, mean, s_dev, ampl); + return ODF_hard(theta, mean, s_dev, ampl) + ODF_hard(theta - simcoon::pi, mean, s_dev, ampl) + ODF_hard(theta + simcoon::pi, mean, s_dev, ampl); break; } case 3: { - return Gaussian(theta, mean, s_dev, ampl) + Gaussian(theta - sim_pi, mean, s_dev, ampl) + Gaussian(theta + sim_pi, mean, s_dev, ampl); + return Gaussian(theta, mean, s_dev, ampl) + Gaussian(theta - simcoon::pi, mean, s_dev, ampl) + Gaussian(theta + simcoon::pi, mean, s_dev, ampl); break; } case 4: { - return Lorentzian(theta, mean, width, ampl) + Lorentzian(theta - sim_pi, mean, width, ampl) + Lorentzian(theta + sim_pi, mean, width, ampl); + return Lorentzian(theta, mean, width, ampl) + Lorentzian(theta - simcoon::pi, mean, width, ampl) + Lorentzian(theta + simcoon::pi, mean, width, ampl); break; } case 5: { - return PseudoVoigt(theta, mean, s_dev, width, ampl, params) + PseudoVoigt(theta - sim_pi, mean, s_dev, width, ampl, params) + PseudoVoigt(theta + sim_pi, mean, s_dev, width, ampl, params); + return PseudoVoigt(theta, mean, s_dev, width, ampl, params) + PseudoVoigt(theta - simcoon::pi, mean, s_dev, width, ampl, params) + PseudoVoigt(theta + simcoon::pi, mean, s_dev, width, ampl, params); break; } case 6: { assert(width > 0.); double inv_width = 1./width; - return Pearson7(theta, mean, inv_width, params) + Pearson7(theta - sim_pi, mean, inv_width, params) + Pearson7(theta + sim_pi, mean, inv_width, params); + return Pearson7(theta, mean, inv_width, params) + Pearson7(theta - simcoon::pi, mean, inv_width, params) + Pearson7(theta + simcoon::pi, mean, inv_width, params); break; } case 7: { diff --git a/src/Continuum_mechanics/Material/read.cpp b/src/Continuum_mechanics/Material/read.cpp index 679531f6b..d5883e58b 100755 --- a/src/Continuum_mechanics/Material/read.cpp +++ b/src/Continuum_mechanics/Material/read.cpp @@ -87,9 +87,9 @@ void read_peak(ODF &odf_rve, const string &path_data, const string &inputfile) { } if(odf_rve.radian == false) { - odf_rve.peaks[i].mean *=(sim_pi/180.); - odf_rve.peaks[i].s_dev *=(sim_pi/180.); - odf_rve.peaks[i].width *=(sim_pi/180.); + odf_rve.peaks[i].mean *=(simcoon::pi/180.); + odf_rve.peaks[i].s_dev *=(simcoon::pi/180.); + odf_rve.peaks[i].width *=(simcoon::pi/180.); } } diff --git a/src/Continuum_mechanics/Micromechanics/schemes.cpp b/src/Continuum_mechanics/Micromechanics/schemes.cpp index d5ff449b7..e91c7c86d 100755 --- a/src/Continuum_mechanics/Micromechanics/schemes.cpp +++ b/src/Continuum_mechanics/Micromechanics/schemes.cpp @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -93,9 +94,14 @@ void Lt_Mori_Tanaka(phase_characteristics &phase, const int &n_matrix) { //Compute the normalization interaction tensir sumT sumT += elli->concentration*elli_multi->T; } - - inv_sumT = inv(sumT); - + + try { + inv_sumT = inv(sumT); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Lt_Mori_Tanaka."); + } + //Compute the strain concentration tensor A for(auto r : phase.sub_phases) { elli_multi = std::dynamic_pointer_cast(r.sptr_multi); @@ -131,8 +137,13 @@ void Lt_Mori_Tanaka_iso(phase_characteristics &phase, const int &n_matrix) { //Compute the normalization interaction tensir sumT sumT += elli->concentration*elli_multi->T; } - - inv_sumT = inv(sumT); + + try { + inv_sumT = inv(sumT); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Lt_Mori_Tanaka_iso."); + } //Compute the strain concentration tensor A for(auto r : phase.sub_phases) { @@ -169,7 +180,12 @@ void DE_Mori_Tanaka(phase_characteristics &phase, const int &n_matrix) { sumT += elli->concentration*elli_multi->T; } - inv_sumT = inv(sumT); + try { + inv_sumT = inv(sumT); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside DE_Mori_Tanaka."); + } //Compute the strain concentration tensor A for(auto r : phase.sub_phases) { @@ -209,7 +225,12 @@ void DE_Mori_Tanaka_iso(phase_characteristics &phase, const int &n_matrix) { sumT += elli->concentration*elli_multi->T; } - inv_sumT = inv(sumT); + try { + inv_sumT = inv(sumT); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside DE_Mori_Tanaka_iso."); + } //Compute the strain concentration tensor A for(auto r : phase.sub_phases) { @@ -382,20 +403,50 @@ void dE_Periodic_Layer(phase_characteristics &phase, const int &nbiter) { lay_multi->sigma_hat(1) = sigma_local(3); lay_multi->sigma_hat(2) = sigma_local(4); } - - for(auto r : phase.sub_phases) { - lay_multi = std::dynamic_pointer_cast(r.sptr_multi); - lay = std::dynamic_pointer_cast(r.sptr_shape); - sumDnn += lay->concentration*inv(lay_multi->Dnn); - sumcDsig += lay->concentration*inv(lay_multi->Dnn)*lay_multi->sigma_hat; - } - vec m = inv(sumDnn)*sumcDsig; - + + try { + for(auto r : phase.sub_phases) { + lay_multi = std::dynamic_pointer_cast(r.sptr_multi); + lay = std::dynamic_pointer_cast(r.sptr_shape); + sumDnn += lay->concentration*inv(lay_multi->Dnn); + sumcDsig += lay->concentration*inv(lay_multi->Dnn)*lay_multi->sigma_hat; + } + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside dE_Periodic_Layer."); + } + + vec m; + try { + m = inv(sumDnn)*sumcDsig; + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside dE_Periodic_Layer."); + } + + try { + for(auto r : phase.sub_phases) { + lay_multi = std::dynamic_pointer_cast(r.sptr_multi); + lay = std::dynamic_pointer_cast(r.sptr_shape); + sumDnn += lay->concentration*inv(lay_multi->Dnn); + sumcDsig += lay->concentration*inv(lay_multi->Dnn)*lay_multi->sigma_hat; + } + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside dE_Periodic_Layer."); + } + for(auto r : phase.sub_phases) { sv_r = std::dynamic_pointer_cast(r.sptr_sv_global); lay_multi = std::dynamic_pointer_cast(r.sptr_multi); - lay_multi->dzdx1 = inv(lay_multi->Dnn)*(m-lay_multi->sigma_hat); + + try { + lay_multi->dzdx1 = inv(lay_multi->Dnn)*(m-lay_multi->sigma_hat); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside dE_Periodic_Layer."); + } lay = std::dynamic_pointer_cast(r.sptr_shape); mat dEtot_local = zeros(6); @@ -449,20 +500,40 @@ void Lt_Periodic_Layer(phase_characteristics &phase) { mat sumDnn = zeros(3,3); mat sumDnt = zeros(3,3); - for(auto r : phase.sub_phases) { - lay_multi = std::dynamic_pointer_cast(r.sptr_multi); - lay = std::dynamic_pointer_cast(r.sptr_shape); - sumDnn += lay->concentration*inv(lay_multi->Dnn); - sumDnt += lay->concentration*inv(lay_multi->Dnn)*lay_multi->Dnt; - } - mat m_n = inv(sumDnn); + + try { + for(auto r : phase.sub_phases) { + lay_multi = std::dynamic_pointer_cast(r.sptr_multi); + lay = std::dynamic_pointer_cast(r.sptr_shape); + sumDnn += lay->concentration*inv(lay_multi->Dnn); + sumDnt += lay->concentration*inv(lay_multi->Dnn)*lay_multi->Dnt; + } + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Lt_Periodic_Layer."); + } + + mat m_n; + try { + m_n = inv(sumDnn); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Lt_Periodic_Layer."); + } + mat m_t = m_n*sumDnt; for(auto r : phase.sub_phases) { lay_multi = std::dynamic_pointer_cast(r.sptr_multi); - lay = std::dynamic_pointer_cast(r.sptr_shape); - lay_multi->dXn = inv(lay_multi->Dnn)*(m_n-lay_multi->Dnn); - lay_multi->dXt = inv(lay_multi->Dnn)*(m_t-lay_multi->Dnt); + lay = std::dynamic_pointer_cast(r.sptr_shape); + + try { + lay_multi->dXn = inv(lay_multi->Dnn)*(m_n-lay_multi->Dnn); + lay_multi->dXt = inv(lay_multi->Dnn)*(m_t-lay_multi->Dnt); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside Lt_Periodic_Layer."); + } A_loc = eye(6,6); diff --git a/src/Continuum_mechanics/Umat/Finite/Neo_hookean_comp.cpp b/src/Continuum_mechanics/Umat/Finite/Neo_hookean_comp.cpp index 8e2472e09..6b8a78f1a 100755 --- a/src/Continuum_mechanics/Umat/Finite/Neo_hookean_comp.cpp +++ b/src/Continuum_mechanics/Umat/Finite/Neo_hookean_comp.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -84,11 +85,23 @@ void umat_neo_hookean_comp(const vec &Etot, const vec &DEtot, const mat &F0, con //Invariants of C double I1 = trace(C); // pow(lambda_alpha(2),2.) + pow(lambda_alpha(1),2.) + pow(lambda_alpha(0),2.); //ascending order - double J = det(F1); //lambda(2)*lambda(1)*lambda(0) + double J; + try { + J = det(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside umat_neo_hookean_comp."); + } double W = mu/2.0*(I1-3.) - mu*log(J) + lambda/2.0*pow(log(J),2.); - mat invC = inv(C); + mat invC; + try { + invC = inv(C); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside umat_neo_hookean_comp."); + } mat I = eye(3,3); //Compute the PKII stress and then the Cauchy stress diff --git a/src/Continuum_mechanics/Umat/Finite/Neo_hookean_incomp.cpp b/src/Continuum_mechanics/Umat/Finite/Neo_hookean_incomp.cpp index ac703c234..26105dc73 100755 --- a/src/Continuum_mechanics/Umat/Finite/Neo_hookean_incomp.cpp +++ b/src/Continuum_mechanics/Umat/Finite/Neo_hookean_incomp.cpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -83,12 +84,24 @@ void umat_neo_hookean_incomp(const vec &etot, const vec &Detot, const mat &F0, c //Invariants of C double I1 = trace(C); // pow(lambda_alpha(2),2.) + pow(lambda_alpha(1),2.) + pow(lambda_alpha(0),2.); //ascending order - double J = det(F1); //lambda(2)*lambda(1)*lambda(0) + double J; + try { + J = det(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside umat_neo_hookean_incomp."); + } double I1_bar = pow(J,-2./3.)*I1; double W = C_10*(I1_bar-3.) + (1./D_1)*pow(J-1.,2.); - mat invC = inv(C); + mat invC; + try { + invC = inv(C); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside umat_neo_hookean_incomp."); + } mat I = eye(3,3); //Compute the PKII stress and then the Cauchy stress diff --git a/src/Continuum_mechanics/Umat/Finite/generic_hyper_invariants.cpp b/src/Continuum_mechanics/Umat/Finite/generic_hyper_invariants.cpp index 263cebb0a..579e42955 100755 --- a/src/Continuum_mechanics/Umat/Finite/generic_hyper_invariants.cpp +++ b/src/Continuum_mechanics/Umat/Finite/generic_hyper_invariants.cpp @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -73,7 +74,13 @@ void umat_generic_hyper_invariants(const std::string &umat_name, const vec &etot double dUdJ = 0.; double dU2dJ2 = 0.; - double J = det(F1); + double J; + try { + J = det(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside umat_generic_hyper_invariants."); + } vec I_bar = isochoric_invariants(b, J); std::map list_potentials; @@ -134,7 +141,7 @@ void umat_generic_hyper_invariants(const std::string &umat_name, const vec &etot double c_2 = props(1); double kappa = props(2); dWdI_1_bar = c_1; - if(fabs(I_bar(1)) > sim_iota) { + if(fabs(I_bar(1)) > simcoon::iota) { dWdI_2_bar = c_2/I_bar(1); dW2dI_22_bar = -1.*c_2/pow(I_bar(1),2.); } diff --git a/src/Continuum_mechanics/Umat/Finite/generic_hyper_pstretch.cpp b/src/Continuum_mechanics/Umat/Finite/generic_hyper_pstretch.cpp index b45483756..9e2c77885 100755 --- a/src/Continuum_mechanics/Umat/Finite/generic_hyper_pstretch.cpp +++ b/src/Continuum_mechanics/Umat/Finite/generic_hyper_pstretch.cpp @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -66,7 +67,13 @@ void umat_generic_hyper_pstretch(const std::string &umat_name, const vec &etot, double dUdJ = 0.; double dU2dJ2 = 0.; - double J = det(F1); + double J; + try { + J = det(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside umat_generic_hyper_pstretch."); + } vec dWdlambda_bar = zeros(3); mat dW2dlambda_bar2 = eye(3,3); diff --git a/src/Continuum_mechanics/Umat/Finite/mooney_rivlin.cpp b/src/Continuum_mechanics/Umat/Finite/mooney_rivlin.cpp index f22455ae1..d3a907727 100755 --- a/src/Continuum_mechanics/Umat/Finite/mooney_rivlin.cpp +++ b/src/Continuum_mechanics/Umat/Finite/mooney_rivlin.cpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -83,12 +84,24 @@ void umat_mooney_rivlin(const vec &Etot, const vec &DEtot, const mat &F0, const //Invariants of C double I1 = trace(C); // pow(lambda_alpha(2),2.) + pow(lambda_alpha(1),2.) + pow(lambda_alpha(0),2.); //ascending order - double J = det(F1); //lambda(2)*lambda(1)*lambda(0) + double J; + try { + J = det(F1); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside umat_mooney_rivlin."); + } double I1_bar = pow(J,-2./3.)*I1; double W = C_10*(I1_bar-3.) + (1./D_1)*pow(J-1.,2.); - mat invC = inv(C); + mat invC; + try { + invC = inv(C); + } catch (const std::runtime_error &e) { + cerr << "Error in inv: " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside umat_mooney_rivlin."); + } mat I = eye(3,3); //Compute the PKII stress and then the Cauchy stress diff --git a/src/Continuum_mechanics/Umat/Mechanical/Combined/Prony_Nfast_Plastic.cpp b/src/Continuum_mechanics/Umat/Mechanical/Combined/Prony_Nfast_Plastic.cpp index 68dbca7a3..30838dc79 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Combined/Prony_Nfast_Plastic.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Combined/Prony_Nfast_Plastic.cpp @@ -138,7 +138,7 @@ void umat_prony_Nfast_plastic(const vec &Etot, const vec &DEtot, vec &sigma, mat double Hp=0.; double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -231,7 +231,7 @@ void umat_prony_Nfast_plastic(const vec &Etot, const vec &DEtot, vec &sigma, mat } //PLastic part - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -254,7 +254,7 @@ void umat_prony_Nfast_plastic(const vec &Etot, const vec &DEtot, vec &sigma, mat Lambdav[i] = eta_norm_strain(flow_visco[i]); dPhi_idv_temp[i] = invH_i[i]*(eta_norm_strain(flow_visco[i])%Ir05()); //Dimension of strain (The flow is of stress type here) - if (DTime > sim_iota) { + if (DTime > simcoon::iota) { Phi(i+1) = norm_strain(flow_visco[i]) - Ds_j(i+1)/DTime; dPhidv[i] = -1.*sum((dPhi_idv_temp[i])%(L_i[i]*Lambdav[i]))-1./DTime; } @@ -319,7 +319,7 @@ void umat_prony_Nfast_plastic(const vec &Etot, const vec &DEtot, vec &sigma, mat for (int i=0; i sim_iota) + if(Ds_j(i) > simcoon::iota) op(i) = 1.; for (int j = 0; j sim_limit) + if (p_ts > simcoon::limit) Hp_ts = beta_ts*pow(p_ts, alpha_ts); else - Hp_ts = sim_iota; + Hp_ts = simcoon::iota; //effective stress sigma_eff = el_pred(L,Eel,ndi); @@ -242,10 +242,10 @@ void umat_damage_LLD_0(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, m //Compute the explicit flow direction Lambdap_ts = eta_stress(sigma_eff_ts); - if (p_ts > sim_limit) + if (p_ts > simcoon::limit) dPhip_tsdp_ts = -1.*alpha_ts*beta_ts*pow(p_ts, alpha_ts-1.); else - dPhip_tsdp_ts = sim_iota; + dPhip_tsdp_ts = simcoon::iota; dPhi_p_tsd_sigma = Theta_ts*eta_stress(sigma_eff_ts); //Here as well @@ -319,22 +319,22 @@ void umat_damage_LLD_0(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, m Eel = Etot + DEtot - alpha*(T+DT-Tinit) - EP; sigma = el_pred(L_tilde, Eel, ndi); - if (fabs(sigma(1)) < sim_iota ) + if (fabs(sigma(1)) < simcoon::iota ) Yd_22 = 0.; else Yd_22 = 0.5*(pow(Macaulay_p(sigma(1)),2.)/(E2_0*pow(1.-d_22,2.))); - if (fabs(sigma(2)) < sim_iota ) + if (fabs(sigma(2)) < simcoon::iota ) Yd_33 = 0.; else Yd_33 = 0.5*(pow(Macaulay_p(sigma(2)),2.)/(E2_0*pow(1.-d_22,2.))); - if (fabs(sigma(3)) < sim_iota ) + if (fabs(sigma(3)) < simcoon::iota ) Yd_12 = 0.; else Yd_12 = 0.5*(pow(sigma(3),2.)/(G12_0*pow(1.-d_12,2.))); - if (fabs(sigma(4)) < sim_iota ) + if (fabs(sigma(4)) < simcoon::iota ) Yd_13 = 0.; else Yd_13 = 0.5*(pow(sigma(4),2.)/(G12_0*pow(1.-d_12,2.))); @@ -347,29 +347,29 @@ void umat_damage_LLD_0(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, m lambda_22 = lagrange_pow_1(d_22, c_lambda, p0_lambda, n_lambda, alpha_lambda); //Preliminaries to compute damage the derivatives of Phi - if (fabs(sigma(1)) < sim_iota ) + if (fabs(sigma(1)) < simcoon::iota ) dYd_22dd = 0.; else dYd_22dd = -0.25*(pow(Macaulay_p(sigma(1)),2.)/(E2_0*pow(1-d_22,3.))); - if (fabs(sigma(2)) < sim_iota ) + if (fabs(sigma(2)) < simcoon::iota ) dYd_33dd = 0.; else dYd_33dd = -0.25*(pow(Macaulay_p(sigma(2)),2.)/(E2_0*pow(1-d_22,3.))); - if (fabs(sigma(3)) < sim_iota ) + if (fabs(sigma(3)) < simcoon::iota ) dYd_12dd = 0.; else dYd_12dd = -0.25*(pow(sigma(3),2.)/(G12_0*pow(1-d_12,3.))); - if (fabs(sigma(4)) < sim_iota ) + if (fabs(sigma(4)) < simcoon::iota ) dYd_13dd = 0.; else dYd_13dd = -0.25*(pow(sigma(4),2.)/(G12_0*pow(1-d_12,3.))); - if (fabs(Yd_12 + Yd_13 + b*(Yd_22 + Yd_33)) < sim_iota ) { + if (fabs(Yd_12 + Yd_13 + b*(Yd_22 + Yd_33)) < simcoon::iota ) { dY_tsdd_22 = 0.; dY_tsdd_12 = 0.; } @@ -390,7 +390,7 @@ void umat_damage_LLD_0(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, m dYts_d_Y12 = 0.; dYts_d_Y13 = 0.; - if (Y_ts > sim_limit) { + if (Y_ts > simcoon::limit) { dYts_d_Y22 = b/(2.*Y_ts); dYts_d_Y33 = b/(2.*Y_ts); dYts_d_Y12 = 1./(2.*Y_ts); @@ -494,7 +494,7 @@ void umat_damage_LLD_0(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, m dYts_d_Y12 = 0.; dYts_d_Y13 = 0.; - if (Y_ts > sim_limit) { + if (Y_ts > simcoon::limit) { dYts_d_Y22 = b/(2.*Y_ts); dYts_d_Y33 = b/(2.*Y_ts); dYts_d_Y12 = 1./(2.*Y_ts); @@ -588,11 +588,11 @@ void umat_damage_LLD_0(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, m vec op = zeros(3); mat delta = eye(3,3); - if(Dd(0) > sim_iota) + if(Dd(0) > simcoon::iota) op(0) = 1.; - if(Dd(1) > sim_iota) + if(Dd(1) > simcoon::iota) op(1) = 1.; - if(Dp(0) > sim_iota) + if(Dp(0) > simcoon::iota) op(2) = 1.; diff --git a/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.cpp b/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.cpp index bb6bac4e4..a0e75f3cd 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.cpp @@ -97,7 +97,7 @@ void umat_damage_weibull(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, double A_Damage = 0.; double Rv = 0.; - if (Mises_stress(sigma) > sim_iota) { + if (Mises_stress(sigma) > simcoon::iota) { Rv = (2./3.)*(1.+nu)+3.*(1.-2.*nu)*pow(tr(sigma)/(3.*Mises_stress(sigma)),2.); A_Damage = pow(Mises_stress(sigma),2.)*Rv/(2.*E*pow(1.-Damage,2.)); } @@ -200,7 +200,7 @@ void umat_damage_weibull(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } @@ -226,7 +226,7 @@ void umat_damage_weibull(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, Lt = L_tilde - (kappa_j[0]*P_epsilon[0].t()); //Lt = L_tilde; - if (Mises_stress(sigma) > sim_iota) { + if (Mises_stress(sigma) > simcoon::iota) { Rv = (2./3.)*(1.+nu)+3.*(1.-2.*nu)*pow(tr(sigma)/(3.*Mises_stress(sigma)),2.); A_Damage = pow(Mises_stress(sigma),2.)*Rv/(2.*E*pow(1.-Damage,2.)); } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.cpp index c9ba568a9..df470ebfd 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.cpp @@ -182,7 +182,7 @@ void umat_ani_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &L //Additional parameters and variables double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = b*(Q-Hp); } else { @@ -239,7 +239,7 @@ void umat_ani_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &L for (compteur = 0; ((compteur < maxiter_umat) && (error > precision_umat)); compteur++) { p = s_j(0); - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = b*(Q-Hp); } else { @@ -297,7 +297,7 @@ void umat_ani_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &L mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.cpp index 28452ebd0..35ed54fae 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.cpp @@ -180,7 +180,7 @@ void umat_dfa_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &L //Additional parameters and variables double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = b*(Q-Hp); } else { @@ -237,7 +237,7 @@ void umat_dfa_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &L for (compteur = 0; ((compteur < maxiter_umat) && (error > precision_umat)); compteur++) { p = s_j(0); - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = b*(Q-Hp); } else { @@ -295,7 +295,7 @@ void umat_dfa_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &L mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.cpp index d5a5b32f5..7b431ecd6 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.cpp @@ -222,7 +222,7 @@ void umat_generic_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, ma //Additional parameters and variables double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { for (int i=0; i precision_umat)); compteur++) { p = s_j(0); - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = 0.; for (int i=0; i sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.cpp index e15d0c295..4d2222c1b 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.cpp @@ -179,7 +179,7 @@ void umat_hill_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat & //Additional parameters and variables double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = b*(Q-Hp); } else { @@ -236,7 +236,7 @@ void umat_hill_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat & for (compteur = 0; ((compteur < maxiter_umat) && (error > precision_umat)); compteur++) { p = s_j(0); - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = b*(Q-Hp); } else { @@ -294,7 +294,7 @@ void umat_hill_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat & mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.cpp index 3865fac2f..e99b1e6ea 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.cpp @@ -147,7 +147,7 @@ void umat_plasticity_hill_isoh_CCP(const vec &Etot, const vec &DEtot, vec &sigma double Hp=0.; double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -194,7 +194,7 @@ void umat_plasticity_hill_isoh_CCP(const vec &Etot, const vec &DEtot, vec &sigma for (compteur = 0; ((compteur < maxiter_umat) && (error > precision_umat)); compteur++) { p = s_j(0); - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -238,7 +238,7 @@ void umat_plasticity_hill_isoh_CCP(const vec &Etot, const vec &DEtot, vec &sigma mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.cpp index 998de2678..6f0095075 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.cpp @@ -145,7 +145,7 @@ void umat_plasticity_hill_isoh_CCP_N(const vec &Etot, const vec &DEtot, vec &sig vec dHpdp = zeros(N_plas); for (int i=0; i sim_iota) { + if (p[i] > simcoon::iota) { dHpdp[i] = m[i]*k[i]*pow(p[i], m[i]-1.0); Hp[i] = k[i]*pow(p[i], m[i]); } @@ -195,7 +195,7 @@ void umat_plasticity_hill_isoh_CCP_N(const vec &Etot, const vec &DEtot, vec &sig p = s_j; for (int i=0; i sim_iota) { + if (p[i] > simcoon::iota) { dHpdp[i] = m[i]*k[i]*pow(p[i], m[i]-1.0); Hp[i] = k[i]*pow(p[i], m[i]); } @@ -248,7 +248,7 @@ void umat_plasticity_hill_isoh_CCP_N(const vec &Etot, const vec &DEtot, vec &sig mat delta = eye(N_plas,N_plas); for (int i=0; i sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.cpp index 2f61984a6..d46fb9709 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.cpp @@ -168,7 +168,7 @@ void umat_plasticity_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, //Additional parameters and variables double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = b*(Q-Hp); } else { @@ -224,7 +224,7 @@ void umat_plasticity_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, for (compteur = 0; ((compteur < maxiter_umat) && (error > precision_umat)); compteur++) { p = s_j(0); - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = b*(Q-Hp); } else { @@ -282,7 +282,7 @@ void umat_plasticity_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.cpp index fa4238780..fa96207dc 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.cpp @@ -138,7 +138,7 @@ void umat_plasticity_iso_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat double Hp=0.; double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -185,7 +185,7 @@ void umat_plasticity_iso_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat for (compteur = 0; ((compteur < maxiter_umat) && (error > precision_umat)); compteur++) { p = s_j(0); - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -231,7 +231,7 @@ void umat_plasticity_iso_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.cpp index 713660794..9a3974b08 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.cpp @@ -135,7 +135,7 @@ void umat_plasticity_kin_iso_CCP(const vec &Etot, const vec &DEtot, vec &sigma, double Hp=0.; double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -188,7 +188,7 @@ void umat_plasticity_kin_iso_CCP(const vec &Etot, const vec &DEtot, vec &sigma, for (compteur = 0; ((compteur < maxiter_umat) && (error > precision_umat)); compteur++) { p = s_j(0); - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -239,7 +239,7 @@ void umat_plasticity_kin_iso_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.cpp b/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.cpp index 4af00c628..7b0244a4b 100644 --- a/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.cpp @@ -22,7 +22,7 @@ #include #include #include - +#include #include #include #include @@ -169,7 +169,7 @@ void umat_sma_mono(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat & ET = vide; xi = 0.; for(int i=0; i 0.) + try { + if(fabs(det(dPhidxi)) > 0.) dxin = -1.*inv(dPhidxi)*Phi; else { //pnewdt = 0.1; // Phi = zeros(nactive); dxin = zeros(nactive); } + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside umat_sma_mono."); + } dxi = 0.; for(int i=0; i 0.)&&(fabs(xi - xi_start) > 1.E-5)) { - if(fabs(det(dPhidxi)) > 0.) { - Lt = L + ((L*lambda)*(inv(dPhidxi))*(dPhidsigma*L)); - } - else - Lt = L; - - vec eigval = eig_sym(Lt); + //Computation of the tangent modulus ! + try { + if(fabs(det(dPhidxi)) > 0.) { + Lt = L + ((L*lambda)*(inv(dPhidxi))*(dPhidsigma*L)); + } + else + Lt = L; + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside umat_sma_mono."); + } + /* ///Note : To use with Self-consistent micromechanics only! + vec eigval = eig_sym(Lt); while(eigval(0) < 0.) { Lt = 0.99*Lt + 0.01*Lt_eff; eigval = eig_sym(Lt); @@ -414,9 +424,8 @@ void umat_sma_mono(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat & Lt = L; } - vec Dsigma = sigma - sigma_start; - vec eigval = eig_sym(Lt); - +// vec Dsigma = sigma - sigma_start; + /* if (compteur < 20) { tnew_dt = 2.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.cpp b/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.cpp index 541b5c4a7..982055e6d 100644 --- a/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.cpp @@ -24,6 +24,7 @@ #include #include +#include #include #include #include @@ -174,7 +175,7 @@ void umat_sma_mono_cubic(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, ET = vide; xi = 0.; for(int i=0; i 0.) + try { + if(fabs(det(dPhidxi)) > 0.) dxin = -1.*inv(dPhidxi)*Phi; else { //pnewdt = 0.1; // Phi = zeros(nactive); dxin = zeros(nactive); } + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside umat_sma_mono_cubic."); + } dxi = 0.; for(int i=0; i 0.)&&(fabs(xi - xi_start) > 1.E-5)) { - if(fabs(det(dPhidxi)) > 0.) { - Lt = L + ((L*lambda)*(inv(dPhidxi))*(dPhidsigma*L)); - } - else - Lt = L; - - vec eigval = eig_sym(Lt); - + try { + if(fabs(det(dPhidxi)) > 0.) { + Lt = L + ((L*lambda)*(inv(dPhidxi))*(dPhidsigma*L)); + } + else + Lt = L; + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_inv("Error in inv function inside umat_sma_mono_cubic."); + } + ///Note : To use with Self-consistent micromechanics only! -/* while(eigval(0) < 0.) { +/* vec eigval = eig_sym(Lt); + while(eigval(0) < 0.) { Lt = 0.99*Lt + 0.01*Lt_eff; eigval = eig_sym(Lt); } */ diff --git a/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp b/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp index b5274a769..44cb7ebcd 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp @@ -193,7 +193,7 @@ void umat_sma_unified_T(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, vec vide = zeros(6); sigma = zeros(6); ET = zeros(6); - xiF = sim_limit; + xiF = simcoon::limit; xiR = 0.; xi = xiF; @@ -609,7 +609,7 @@ void umat_sma_unified_T(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat delta = eye(2,2); for (int i=0; i<2; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.cpp b/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.cpp index 488ced93c..c1643a341 100644 --- a/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.cpp @@ -184,7 +184,7 @@ void umat_prony_Nfast(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, co Lambdav[i] = eta_norm_strain(flow_visco[i]); dPhi_idv_temp[i] = invH_i[i]*(eta_norm_strain(flow_visco[i])%Ir05()); //Dimension of strain (The flow is of stress type here) - if (DTime > sim_iota) { + if (DTime > simcoon::iota) { Phi(i) = norm_strain(flow_visco[i]) - Ds_j(i)/DTime; dPhidv[i] = -1.*sum((dPhi_idv_temp[i])%(L_i[i]*Lambdav[i]))-1./DTime; } @@ -241,7 +241,7 @@ void umat_prony_Nfast(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, co for (int i=0; i sim_iota) + if(Ds_j(i) > simcoon::iota) op(i) = 1.; for (int j = 0; j sim_iota) { + if (DTime > simcoon::iota) { Phi(i) = norm_strain(flow_visco[i]) - Ds_j(i)/DTime; dPhidv[i] = -1.*sum((dPhi_idsigma[i])%(L_i[i]*Lambdav[i]))-1./DTime; } @@ -243,7 +243,7 @@ void umat_zener_Nfast(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, co for (int i=0; i sim_iota) + if(Ds_j(i) > simcoon::iota) op(i) = 1.; for (int j = 0; j sim_iota) { + if (DTime > simcoon::iota) { Phi(0) = norm_strain(flow_V1) - Ds_j(0)/DTime; dPhidv = -1.*sum((dPhidsigma)%(L1*Lambdav))-1./DTime; } @@ -169,7 +169,7 @@ void umat_zener_fast(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, con mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Thermomechanical/Plasticity/plastic_isotropic_ccp.cpp b/src/Continuum_mechanics/Umat/Thermomechanical/Plasticity/plastic_isotropic_ccp.cpp index 0ee9cf82e..58bf2273d 100755 --- a/src/Continuum_mechanics/Umat/Thermomechanical/Plasticity/plastic_isotropic_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Thermomechanical/Plasticity/plastic_isotropic_ccp.cpp @@ -148,7 +148,7 @@ void umat_plasticity_iso_CCP_T(const vec &Etot, const vec &DEtot, vec &sigma, do double Hp=0.; double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -195,7 +195,7 @@ void umat_plasticity_iso_CCP_T(const vec &Etot, const vec &DEtot, vec &sigma, do for (compteur = 0; ((compteur < maxiter_umat) && (error > precision_umat)); compteur++) { p = s_j(0); - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -239,7 +239,7 @@ void umat_plasticity_iso_CCP_T(const vec &Etot, const vec &DEtot, vec &sigma, do mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Thermomechanical/Plasticity/plastic_kin_iso_ccp.cpp b/src/Continuum_mechanics/Umat/Thermomechanical/Plasticity/plastic_kin_iso_ccp.cpp index 2c9feb835..067800c11 100755 --- a/src/Continuum_mechanics/Umat/Thermomechanical/Plasticity/plastic_kin_iso_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Thermomechanical/Plasticity/plastic_kin_iso_ccp.cpp @@ -144,7 +144,7 @@ void umat_plasticity_kin_iso_CCP_T(const vec &Etot, const vec &DEtot, vec &sigma double Hp=0.; double dHpdp=0.; - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -197,7 +197,7 @@ void umat_plasticity_kin_iso_CCP_T(const vec &Etot, const vec &DEtot, vec &sigma for (compteur = 0; ((compteur < maxiter_umat) && (error > precision_umat)); compteur++) { p = s_j(0); - if (p > sim_iota) { + if (p > simcoon::iota) { dHpdp = m*k*pow(p, m-1); Hp = k*pow(p, m); } @@ -246,7 +246,7 @@ void umat_plasticity_kin_iso_CCP_T(const vec &Etot, const vec &DEtot, vec &sigma mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Thermomechanical/SMA/unified_T.cpp b/src/Continuum_mechanics/Umat/Thermomechanical/SMA/unified_T.cpp index 6279d46b6..c2fc5c000 100755 --- a/src/Continuum_mechanics/Umat/Thermomechanical/SMA/unified_T.cpp +++ b/src/Continuum_mechanics/Umat/Thermomechanical/SMA/unified_T.cpp @@ -197,7 +197,7 @@ void umat_sma_unified_T_T(const vec &Etot, const vec &DEtot, vec &sigma, double vec vide = zeros(6); sigma = zeros(6); ET = zeros(6); - xiF = sim_limit; + xiF = simcoon::limit; xiR = 0.; xi = xiF; @@ -601,7 +601,7 @@ void umat_sma_unified_T_T(const vec &Etot, const vec &DEtot, vec &sigma, double mat delta = eye(2,2); for (int i=0; i<2; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/Thermomechanical/Viscoelasticity/Prony_Nfast.cpp b/src/Continuum_mechanics/Umat/Thermomechanical/Viscoelasticity/Prony_Nfast.cpp index 7467493bd..654c4a8cc 100644 --- a/src/Continuum_mechanics/Umat/Thermomechanical/Viscoelasticity/Prony_Nfast.cpp +++ b/src/Continuum_mechanics/Umat/Thermomechanical/Viscoelasticity/Prony_Nfast.cpp @@ -201,7 +201,7 @@ void umat_prony_Nfast_T(const vec &Etot, const vec &DEtot, vec &sigma, double &r dPhi_idv_temp[i] = invH_i[i]*(eta_norm_strain(flow_visco[i])%Ir05()); //Dimension of strain (The flow is of stress type here) kappa_j[i] = L_i[i]*Lambdav[i]; - if (DTime > sim_iota) { + if (DTime > simcoon::iota) { Phi(i) = norm_strain(flow_visco[i]) - Ds_j(i)/DTime; dPhidv[i] = -1.*sum((dPhi_idv_temp[i])%(L_i[i]*Lambdav[i]))-1./DTime; } @@ -271,7 +271,7 @@ void umat_prony_Nfast_T(const vec &Etot, const vec &DEtot, vec &sigma, double &r for (int i=0; i sim_iota) + if(Ds_j(i) > simcoon::iota) op(i) = 1.; for (int j = 0; j sim_iota) { + if (DTime > simcoon::iota) { Phi(i) = norm_strain(flow_visco[i]) - Ds_j(i)/DTime; dPhidv[i] = -1.*sum((dPhi_idsigma[i])%(L_i[i]*Lambdav[i]))-1./DTime; } @@ -276,7 +276,7 @@ void umat_zener_Nfast_T(const vec &Etot, const vec &DEtot, vec &sigma, double &r for (int i=0; i sim_iota) + if(Ds_j(i) > simcoon::iota) op(i) = 1.; for (int j = 0; j sim_iota) { + if (DTime > simcoon::iota) { Phi(0) = norm_stress(sigma_tildeV1) - Ds_j(0)/DTime; dPhidv = -1.*sum((dPhidsigma%Ir2())%(L1*Lambdav))-1./DTime; } @@ -200,7 +200,7 @@ void umat_zener_fast_T(const vec &Etot, const vec &DEtot, vec &sigma, double &r, mat delta = eye(1,1); for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) + if(Ds_j[i] > simcoon::iota) op(i) = 1.; } diff --git a/src/Continuum_mechanics/Umat/umat_smart.cpp b/src/Continuum_mechanics/Umat/umat_smart.cpp index 1a5974c0b..f23357bb5 100644 --- a/src/Continuum_mechanics/Umat/umat_smart.cpp +++ b/src/Continuum_mechanics/Umat/umat_smart.cpp @@ -456,7 +456,7 @@ void run_umat_T(phase_characteristics &rve, const mat &DR,const double &Time,con tnew_dt = 1.; - if (Time > sim_limit) { + if (Time > simcoon::limit) { start = false; } switch (control_type) { @@ -479,7 +479,7 @@ void run_umat_M(phase_characteristics &rve, const mat &DR, const double &Time, c { tnew_dt = 1.; - if (Time > sim_limit) { + if (Time > simcoon::limit) { start = false; } switch (control_type) { diff --git a/src/Simulation/Identification/identification.cpp b/src/Simulation/Identification/identification.cpp index 1e742e954..6123a4af0 100755 --- a/src/Simulation/Identification/identification.cpp +++ b/src/Simulation/Identification/identification.cpp @@ -278,7 +278,7 @@ void run_identification(const std::string &simul_type, const int &n_param, const find_best(gen[g], gboys[g], gen[g-1], gboys[g-1], gensons, maxpop, n_param, id0); write_results(result, outputfile, gen[g], g, maxpop, n_param); - if(fabs(costnm1 - gen[g].pop[0].cout) < sim_iota) { + if(fabs(costnm1 - gen[g].pop[0].cout) < simcoon::iota) { compt_des++; } else{ diff --git a/src/Simulation/Identification/optimize.cpp b/src/Simulation/Identification/optimize.cpp index 12dfc1153..5b102b6e2 100755 --- a/src/Simulation/Identification/optimize.cpp +++ b/src/Simulation/Identification/optimize.cpp @@ -81,7 +81,7 @@ Col checkS(const mat &S) { for(int i=0; i sim_iota) { + if (W(z) > simcoon::iota) { Cout += pow((vexp(z)-vnum(z)), 2.)*W(z); } } diff --git a/src/Simulation/Maths/num_solve.cpp b/src/Simulation/Maths/num_solve.cpp index d20ae0a62..25a8b92d9 100755 --- a/src/Simulation/Maths/num_solve.cpp +++ b/src/Simulation/Maths/num_solve.cpp @@ -19,6 +19,7 @@ #include #include #include +#include using namespace std; using namespace arma; @@ -37,12 +38,17 @@ void Newton_Raphon(const vec &Phi, const vec &Y_crit, const mat &denom, vec &Dp, assert(fabs(Y_crit(i)) > 0.); } - if (fabs(det(denom)) > sim_limit) { - dp = -1.*solve(denom,Phi); - } - else { - dp = zeros(n); - } + try { + if (fabs(det(denom)) > simcoon::limit) { + dp = -1.*solve(denom,Phi); + } + else { + dp = zeros(n); + } + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside Newton_Raphon."); + } //New update of the vector of the variable increment throughout the step Dp = Dp + dp; @@ -100,12 +106,17 @@ void Fischer_Burmeister(const vec &Phi, const vec &Y_crit, const mat &denom, vec } - if (fabs(det(denomFB)) > sim_limit) { - dp = -1.*solve(denomFB,FB); - } - else { - dp = zeros(n); - } + try { + if (fabs(det(denomFB)) > simcoon::limit) { + dp = -1.*solve(denomFB,FB); + } + else { + dp = zeros(n); + } + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside Fischer_Burmeister."); + } //New update of the vector of the variable increment throughout the step Dp = Dp + dp; @@ -218,13 +229,18 @@ void Fischer_Burmeister_limits(const vec &Phi_p, const vec &Phi_l, const vec &Y_ } } - if (fabs(det(denomFB)) > sim_limit) { - dp = -1.*solve(denomFB,FB); - } - else { - dp = zeros(n); - } - + try { + if (fabs(det(denomFB)) > simcoon::limit) { + dp = -1.*solve(denomFB,FB); + } + else { + dp = zeros(n); + } + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside Fischer_Burmeister_limits."); + } + //New update of the vector of the variable increment throughout the step Dp = Dp + dp; error = 0.; @@ -291,12 +307,17 @@ void Fischer_Burmeister_m(const vec &Phi, const vec &Y_crit, const mat &denom, v } - if (fabs(det(denomFB)) > sim_limit) { - dp = -1.*solve(denomFB, FB); - } - else { - dp = zeros(n); - } + try { + if (fabs(det(denomFB)) > simcoon::limit) { + dp = -1.*solve(denomFB, FB); + } + else { + dp = zeros(n); + } + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside Fischer_Burmeister_m."); + } //New update of the transformation/orientation multipliers Dp = Dp + dp; @@ -421,12 +442,17 @@ void Fischer_Burmeister_m_limits(const vec &Phi_p, const vec &Phi_l, const vec & } } - if (fabs(det(denomFB)) > sim_limit) { - dp = -1.*solve(denomFB, FB); - } - else { - dp = zeros(n); - } + try { + if (fabs(det(denomFB)) > simcoon::limit) { + dp = -1.*solve(denomFB, FB); + } + else { + dp = zeros(n); + } + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside Fischer_Burmeister_m_limits."); + } //New update of the transformation/orientation multipliers Dp = Dp + dp; diff --git a/src/Simulation/Maths/rotation.cpp b/src/Simulation/Maths/rotation.cpp index df1e1bbb6..61f0ee3e7 100755 --- a/src/Simulation/Maths/rotation.cpp +++ b/src/Simulation/Maths/rotation.cpp @@ -660,13 +660,13 @@ vec rotate_strain(const vec &V, const mat &DR, const bool &active) { mat rotate_l2g_strain(const vec &E, const double &psi, const double &theta, const double &phi) { mat E_temp = E; - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { E_temp = rotate_strain(E_temp, -phi, axis_phi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { E_temp = rotate_strain(E_temp, -theta, axis_theta, false); } - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { E_temp = rotate_strain(E_temp, -psi, axis_psi, false); } @@ -677,13 +677,13 @@ mat rotate_l2g_strain(const vec &E, const double &psi, const double &theta, cons mat rotate_g2l_strain(const vec &E, const double &psi, const double &theta, const double &phi) { mat E_temp = E; - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { E_temp = rotate_strain(E_temp, psi, axis_psi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { E_temp = rotate_strain(E_temp, theta, axis_theta, false); } - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { E_temp = rotate_strain(E_temp, phi, axis_phi, false); } @@ -694,13 +694,13 @@ mat rotate_g2l_strain(const vec &E, const double &psi, const double &theta, cons mat rotate_l2g_stress(const vec &S, const double &psi, const double &theta, const double &phi) { mat S_temp = S; - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { S_temp = rotate_stress(S_temp, -phi, axis_phi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { S_temp = rotate_stress(S_temp, -theta, axis_theta, false); } - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { S_temp = rotate_stress(S_temp, -psi, axis_psi, false); } @@ -711,13 +711,13 @@ mat rotate_l2g_stress(const vec &S, const double &psi, const double &theta, cons mat rotate_g2l_stress(const vec &S, const double &psi, const double &theta, const double &phi) { mat S_temp = S; - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { S_temp = rotate_stress(S_temp, psi, axis_psi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { S_temp = rotate_stress(S_temp, theta, axis_theta, false); } - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { S_temp = rotate_stress(S_temp, phi, axis_phi, false); } @@ -728,13 +728,13 @@ mat rotate_g2l_stress(const vec &S, const double &psi, const double &theta, cons mat rotate_l2g_L(const mat &Lt, const double &psi, const double &theta, const double &phi) { mat Lt_temp = Lt; - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { Lt_temp = rotateL(Lt_temp, -phi, axis_phi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { Lt_temp = rotateL(Lt_temp, -theta, axis_theta, false); } - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { Lt_temp = rotateL(Lt_temp, -psi, axis_psi, false); } @@ -745,13 +745,13 @@ mat rotate_l2g_L(const mat &Lt, const double &psi, const double &theta, const do mat rotate_g2l_L(const mat &Lt, const double &psi, const double &theta, const double &phi) { mat Lt_temp = Lt; - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { Lt_temp = rotateL(Lt_temp, psi, axis_psi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { Lt_temp = rotateL(Lt_temp, theta, axis_theta, false); } - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { Lt_temp = rotateL(Lt_temp, phi, axis_phi, false); } @@ -762,13 +762,13 @@ mat rotate_g2l_L(const mat &Lt, const double &psi, const double &theta, const do mat rotate_l2g_M(const mat &M, const double &psi, const double &theta, const double &phi) { mat M_temp = M; - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { M_temp = rotateM(M_temp, -phi, axis_phi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { M_temp = rotateM(M_temp, -theta, axis_theta, false); } - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { M_temp = rotateM(M_temp, -psi, axis_psi, false); } @@ -779,13 +779,13 @@ mat rotate_l2g_M(const mat &M, const double &psi, const double &theta, const dou mat rotate_g2l_M(const mat &M, const double &psi, const double &theta, const double &phi) { mat M_temp = M; - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { M_temp = rotateM(M_temp, psi, axis_psi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { M_temp = rotateM(M_temp, theta, axis_theta, false); } - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { M_temp = rotateM(M_temp, phi, axis_phi, false); } @@ -796,13 +796,13 @@ mat rotate_g2l_M(const mat &M, const double &psi, const double &theta, const dou mat rotate_l2g_A(const mat &A, const double &psi, const double &theta, const double &phi) { mat A_temp = A; - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { A_temp = rotateA(A_temp, -phi, axis_phi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { A_temp = rotateA(A_temp, -theta, axis_theta, false); } - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { A_temp = rotateA(A_temp, -psi, axis_psi, false); } @@ -813,13 +813,13 @@ mat rotate_l2g_A(const mat &A, const double &psi, const double &theta, const dou mat rotate_g2l_A(const mat &A, const double &psi, const double &theta, const double &phi) { mat A_temp = A; - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { A_temp = rotateA(A_temp, psi, axis_psi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { A_temp = rotateA(A_temp, theta, axis_theta, false); } - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { A_temp = rotateA(A_temp, phi, axis_phi, false); } @@ -830,13 +830,13 @@ mat rotate_g2l_A(const mat &A, const double &psi, const double &theta, const dou mat rotate_l2g_B(const mat &B, const double &psi, const double &theta, const double &phi) { mat B_temp = B; - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { B_temp = rotateB(B_temp, -phi, axis_phi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { B_temp = rotateB(B_temp, -theta, axis_theta, false); } - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { B_temp = rotateB(B_temp, -psi, axis_psi, false); } @@ -847,13 +847,13 @@ mat rotate_l2g_B(const mat &B, const double &psi, const double &theta, const dou mat rotate_g2l_B(const mat &B, const double &psi, const double &theta, const double &phi) { mat B_temp = B; - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { B_temp = rotateB(B_temp, psi, axis_psi, false); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { B_temp = rotateB(B_temp, theta, axis_theta, false); } - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { B_temp = rotateB(B_temp, phi, axis_phi, false); } diff --git a/src/Simulation/Maths/stats.cpp b/src/Simulation/Maths/stats.cpp index aa1d59827..32b8c49a9 100755 --- a/src/Simulation/Maths/stats.cpp +++ b/src/Simulation/Maths/stats.cpp @@ -41,7 +41,7 @@ double normal_distrib(const double &x, const double &mean, const double &dev){ else { double k = 1.0/(1.0 + 0.2316419*x_norm); double k_sum = k*(0.319381530 + k*(-0.356563782 + k*(1.781477937 + k*(-1.821255978 + 1.330274429*k)))); - return (1.0 - (1.0/(sqrt(2*sim_pi)))*exp(-0.5*x_norm*x_norm) * k_sum); + return (1.0 - (1.0/(sqrt(2*simcoon::pi)))*exp(-0.5*x_norm*x_norm) * k_sum); } } @@ -86,7 +86,7 @@ double ODF_sd(const double &Theta, const double &mean, const vec ¶ms) { if (fabs(Theta - mean) < 1.E-6) return alpha1; - else if (fabs((Theta - mean)-0.5*sim_pi) < 1.E-6) + else if (fabs((Theta - mean)-0.5*simcoon::pi) < 1.E-6) return 0.; else return fabs((alpha1*pow(cos(Theta - mean),2.*pow1) + alpha2*pow(cos(Theta - mean),2.*pow2)*pow(sin(Theta - mean),2.*pow2))*cos(Theta - mean)); @@ -106,7 +106,7 @@ double Gaussian(const double &X, const double &mean, const double &std_dev, cons assert(std_dev>0); - return ampl/(std_dev * sqrt(2.*sim_pi)) * exp(-1./2. * pow((X - mean)/std_dev, 2.)); + return ampl/(std_dev * sqrt(2.*simcoon::pi)) * exp(-1./2. * pow((X - mean)/std_dev, 2.)); } ///4. Lorentzian @@ -114,7 +114,7 @@ double Lorentzian(const double &X, const double &mean, const double &width, cons assert(width>0); - return ampl * width/(2.*sim_pi*(pow(X - mean, 2.)+ pow(width / 2., 2.))); + return ampl * width/(2.*simcoon::pi*(pow(X - mean, 2.)+ pow(width / 2., 2.))); } ///5. Pseudo-Voigt @@ -134,7 +134,7 @@ double Pearson7(const double& X, const double& mean, const double &inv_width, ve double max = params(0); double shape = params(1); assert(shape>0); - if (fabs(max) < sim_limit) { + if (fabs(max) < simcoon::limit) { max = 1.; } diff --git a/src/Simulation/Phase/phase_characteristics.cpp b/src/Simulation/Phase/phase_characteristics.cpp index a2d0bdf56..3c1a94da4 100755 --- a/src/Simulation/Phase/phase_characteristics.cpp +++ b/src/Simulation/Phase/phase_characteristics.cpp @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -417,7 +418,13 @@ void phase_characteristics::output(const solver_output &so, const int &kblock, c break; } case 1: { - vec E_biot = t2v_strain(sqrtmat_sympd(2.*v2t_strain(sptr_sv_global->Etot)+eye(3,3)) - eye(3,3)); + vec E_biot; + try { + E_biot = t2v_strain(sqrtmat_sympd(2.*v2t_strain(sptr_sv_global->Etot)+eye(3,3)) - eye(3,3)); + } catch (const std::runtime_error &e) { + cerr << "Error in sqrtmat_sympd : " << e.what() << endl; + throw simcoon::exception_sqrtmat_sympd("Error in sqrtmat_sympd function inside phase_characteristics::output."); + } for (int z=0; zEtot)+eye(3,3)) - eye(3,3)); + vec E_biot; + try { + E_biot = t2v_strain(sqrtmat_sympd(2.*v2t_strain(sptr_sv_local->Etot)+eye(3,3)) - eye(3,3)); + } catch (const std::runtime_error &e) { + cerr << "Error in sqrtmat_sympd : " << e.what() << endl; + throw simcoon::exception_sqrtmat_sympd("Error in sqrtmat_sympd function inside phase_characteristics::output."); + } for (int z=0; z> r.sptr_matprops->props(j); } - r.sptr_matprops->psi_mat*=(sim_pi/180.); - r.sptr_matprops->theta_mat*=(sim_pi/180.); - r.sptr_matprops->phi_mat*=(sim_pi/180.); + r.sptr_matprops->psi_mat*=(simcoon::pi/180.); + r.sptr_matprops->theta_mat*=(simcoon::pi/180.); + r.sptr_matprops->phi_mat*=(simcoon::pi/180.); } paramphases.close(); @@ -210,13 +210,13 @@ void read_layer(phase_characteristics &rve, const string &path_data, const strin paramphases >> r.sptr_matprops->props(j); } - r.sptr_matprops->psi_mat*=(sim_pi/180.); - r.sptr_matprops->theta_mat*=(sim_pi/180.); - r.sptr_matprops->phi_mat*=(sim_pi/180.); + r.sptr_matprops->psi_mat*=(simcoon::pi/180.); + r.sptr_matprops->theta_mat*=(simcoon::pi/180.); + r.sptr_matprops->phi_mat*=(simcoon::pi/180.); - sptr_layer->psi_geom*=(sim_pi/180.); - sptr_layer->theta_geom*=(sim_pi/180.); - sptr_layer->phi_geom*=(sim_pi/180.); + sptr_layer->psi_geom*=(simcoon::pi/180.); + sptr_layer->theta_geom*=(simcoon::pi/180.); + sptr_layer->phi_geom*=(simcoon::pi/180.); } paramphases.close(); @@ -284,13 +284,13 @@ void read_ellipsoid(phase_characteristics &rve, const string &path_data, const s paramphases >> r.sptr_matprops->props(j); } - r.sptr_matprops->psi_mat*=(sim_pi/180.); - r.sptr_matprops->theta_mat*=(sim_pi/180.); - r.sptr_matprops->phi_mat*=(sim_pi/180.); + r.sptr_matprops->psi_mat*=(simcoon::pi/180.); + r.sptr_matprops->theta_mat*=(simcoon::pi/180.); + r.sptr_matprops->phi_mat*=(simcoon::pi/180.); - sptr_ellipsoid->psi_geom*=(sim_pi/180.); - sptr_ellipsoid->theta_geom*=(sim_pi/180.); - sptr_ellipsoid->phi_geom*=(sim_pi/180.); + sptr_ellipsoid->psi_geom*=(simcoon::pi/180.); + sptr_ellipsoid->theta_geom*=(simcoon::pi/180.); + sptr_ellipsoid->phi_geom*=(simcoon::pi/180.); } paramphases.close(); @@ -369,13 +369,13 @@ void read_cylinder(phase_characteristics &rve, const string &path_data, const st paramphases >> r.sptr_matprops->props(j); } - r.sptr_matprops->psi_mat*=(sim_pi/180.); - r.sptr_matprops->theta_mat*=(sim_pi/180.); - r.sptr_matprops->phi_mat*=(sim_pi/180.); + r.sptr_matprops->psi_mat*=(simcoon::pi/180.); + r.sptr_matprops->theta_mat*=(simcoon::pi/180.); + r.sptr_matprops->phi_mat*=(simcoon::pi/180.); - sptr_cylinder->psi_geom*=(sim_pi/180.); - sptr_cylinder->theta_geom*=(sim_pi/180.); - sptr_cylinder->phi_geom*=(sim_pi/180.); + sptr_cylinder->psi_geom*=(simcoon::pi/180.); + sptr_cylinder->theta_geom*=(simcoon::pi/180.); + sptr_cylinder->phi_geom*=(simcoon::pi/180.); } paramphases.close(); diff --git a/src/Simulation/Phase/state_variables.cpp b/src/Simulation/Phase/state_variables.cpp index 90d45b148..38a97c921 100755 --- a/src/Simulation/Phase/state_variables.cpp +++ b/src/Simulation/Phase/state_variables.cpp @@ -498,7 +498,7 @@ state_variables& state_variables::rotate_l2g(const state_variables& sv, const do statev = sv.statev; statev_start = sv.statev_start; - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { Etot = rotate_strain(Etot, -phi, axis_phi); DEtot = rotate_strain(DEtot, -phi, axis_phi); etot = rotate_strain(etot, -phi, axis_phi); @@ -516,7 +516,7 @@ state_variables& state_variables::rotate_l2g(const state_variables& sv, const do R = rotate_mat(R, -phi, axis_phi); DR = rotate_mat(DR, -phi, axis_phi); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { Etot = rotate_strain(Etot, -theta, axis_theta); DEtot = rotate_strain(DEtot, -theta, axis_theta); etot = rotate_strain(etot, -theta, axis_theta); @@ -534,7 +534,7 @@ state_variables& state_variables::rotate_l2g(const state_variables& sv, const do R = rotate_mat(R, -theta, axis_theta); DR = rotate_mat(DR, -theta, axis_theta); } - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { Etot = rotate_strain(Etot, -psi, axis_psi); DEtot = rotate_strain(DEtot, -psi, axis_psi); etot = rotate_strain(etot, -psi, axis_psi); @@ -583,7 +583,7 @@ state_variables& state_variables::rotate_g2l(const state_variables& sv, const do statev = sv.statev; statev_start = sv.statev_start; - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { Etot = rotate_strain(Etot, psi, axis_psi); DEtot = rotate_strain(DEtot, psi, axis_psi); etot = rotate_strain(etot, psi, axis_psi); @@ -601,7 +601,7 @@ state_variables& state_variables::rotate_g2l(const state_variables& sv, const do R = rotate_mat(R, psi, axis_psi); DR = rotate_mat(DR, psi, axis_psi); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { Etot = rotate_strain(Etot, theta, axis_theta); DEtot = rotate_strain(DEtot, theta, axis_theta); etot = rotate_strain(etot, theta, axis_theta); @@ -619,7 +619,7 @@ state_variables& state_variables::rotate_g2l(const state_variables& sv, const do R = rotate_mat(R, theta, axis_theta); DR = rotate_mat(DR, theta, axis_theta); } - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { Etot = rotate_strain(Etot, phi, axis_phi); DEtot = rotate_strain(DEtot, phi, axis_phi); etot = rotate_strain(etot, phi, axis_phi); diff --git a/src/Simulation/Phase/state_variables_M.cpp b/src/Simulation/Phase/state_variables_M.cpp index 4d134f594..3c9f045ee 100755 --- a/src/Simulation/Phase/state_variables_M.cpp +++ b/src/Simulation/Phase/state_variables_M.cpp @@ -249,19 +249,19 @@ state_variables_M& state_variables_M::rotate_l2g(const state_variables_M& sv, co L = sv.L; Lt = sv.Lt; - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, -phi, axis_phi); sigma_in_start = rotate_stress(sigma_in_start, -phi, axis_phi); L = rotateL(L, -phi, axis_phi); Lt = rotateL(Lt, -phi, axis_phi); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, -theta, axis_theta); sigma_in_start = rotate_stress(sigma_in_start, -theta, axis_theta); Lt = rotateL(Lt, -theta, axis_theta); L = rotateL(L, -theta, axis_theta); } - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, -psi, axis_psi); sigma_in_start = rotate_stress(sigma_in_start, -psi, axis_psi); Lt = rotateL(Lt, -psi, axis_psi); @@ -287,19 +287,19 @@ state_variables_M& state_variables_M::rotate_g2l(const state_variables_M& sv, co L = sv.L; Lt = sv.Lt; - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, psi, axis_psi); sigma_in_start = rotate_stress(sigma_in_start, psi, axis_psi); L = rotateL(L, psi, axis_psi); Lt = rotateL(Lt, psi, axis_psi); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, theta, axis_theta); sigma_in_start = rotate_stress(sigma_in_start, theta, axis_theta); L = rotateL(L, theta, axis_theta); Lt = rotateL(Lt, theta, axis_theta); } - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, phi, axis_phi); sigma_in_start = rotate_stress(sigma_in_start, phi, axis_phi); L = rotateL(L, phi, axis_phi); diff --git a/src/Simulation/Phase/state_variables_T.cpp b/src/Simulation/Phase/state_variables_T.cpp index ebd513341..2572aa399 100755 --- a/src/Simulation/Phase/state_variables_T.cpp +++ b/src/Simulation/Phase/state_variables_T.cpp @@ -321,7 +321,7 @@ state_variables_T& state_variables_T::rotate_l2g(const state_variables_T& sv, co Wm_start = sv.Wm_start; Wt_start = sv.Wt_start; - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, -phi, axis_phi); sigma_in_start = rotate_stress(sigma_in_start, -phi, axis_phi); dSdE = rotateL(dSdE, -phi, axis_phi); @@ -329,7 +329,7 @@ state_variables_T& state_variables_T::rotate_l2g(const state_variables_T& sv, co dSdT = rotate_stress(dSdT, -phi, axis_phi); drdE = rotate_strain(drdE, -phi, axis_phi); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, -theta, axis_theta); sigma_in_start = rotate_stress(sigma_in_start, -theta, axis_theta); dSdE = rotateL(dSdE, -theta, axis_theta); @@ -337,7 +337,7 @@ state_variables_T& state_variables_T::rotate_l2g(const state_variables_T& sv, co dSdT = rotate_stress(dSdT, -theta, axis_theta); drdE = rotate_strain(drdE, -theta, axis_theta); } - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, -psi, axis_psi); sigma_in_start = rotate_stress(sigma_in_start, -psi, axis_psi); dSdE = rotateL(dSdE, -psi, axis_psi); @@ -371,7 +371,7 @@ state_variables_T& state_variables_T::rotate_g2l(const state_variables_T& sv, co Wm_start = sv.Wm_start; Wt_start = sv.Wt_start; - if(fabs(psi) > sim_iota) { + if(fabs(psi) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, psi, axis_psi); sigma_in_start = rotate_stress(sigma_in_start, psi, axis_psi); dSdE = rotateL(dSdE, psi, axis_psi); @@ -380,7 +380,7 @@ state_variables_T& state_variables_T::rotate_g2l(const state_variables_T& sv, co drdE = rotate_strain(drdE, psi, axis_psi); } - if(fabs(theta) > sim_iota) { + if(fabs(theta) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, theta, axis_theta); sigma_in_start = rotate_stress(sigma_in_start, theta, axis_theta); dSdE = rotateL(dSdE, theta, axis_theta); @@ -388,7 +388,7 @@ state_variables_T& state_variables_T::rotate_g2l(const state_variables_T& sv, co dSdT = rotate_stress(dSdT, theta, axis_theta); drdE = rotate_strain(drdE, theta, axis_theta); } - if(fabs(phi) > sim_iota) { + if(fabs(phi) > simcoon::iota) { sigma_in = rotate_stress(sigma_in, phi, axis_phi); sigma_in_start = rotate_stress(sigma_in_start, phi, axis_phi); dSdE = rotateL(dSdE, phi, axis_phi); diff --git a/src/Simulation/Phase/write.cpp b/src/Simulation/Phase/write.cpp index e8f79c6db..56b5ba0bb 100755 --- a/src/Simulation/Phase/write.cpp +++ b/src/Simulation/Phase/write.cpp @@ -48,9 +48,9 @@ void write_phase(phase_characteristics &rve, const string &path_data, const stri for(auto r : rve.sub_phases) { - r.sptr_matprops->psi_mat*=(180./sim_pi); - r.sptr_matprops->theta_mat*=(180./sim_pi); - r.sptr_matprops->phi_mat*=(180./sim_pi); + r.sptr_matprops->psi_mat*=(180./simcoon::pi); + r.sptr_matprops->theta_mat*=(180./simcoon::pi); + r.sptr_matprops->phi_mat*=(180./simcoon::pi); paramphases << r.sptr_matprops->number << "\t" << r.sptr_matprops->umat_name << "\t" << r.sptr_matprops->save << "\t" << r.sptr_shape->concentration << "\t" << r.sptr_matprops->psi_mat << "\t" << r.sptr_matprops->theta_mat << "\t" << r.sptr_matprops->phi_mat << "\t" << r.sptr_matprops->nprops << "\t" << r.sptr_sv_global->nstatev; @@ -74,14 +74,14 @@ void write_layer(phase_characteristics &rve, const string &path_data, const stri for(auto r : rve.sub_phases) { - r.sptr_matprops->psi_mat*=(180./sim_pi); - r.sptr_matprops->theta_mat*=(180./sim_pi); - r.sptr_matprops->phi_mat*=(180./sim_pi); + r.sptr_matprops->psi_mat*=(180./simcoon::pi); + r.sptr_matprops->theta_mat*=(180./simcoon::pi); + r.sptr_matprops->phi_mat*=(180./simcoon::pi); auto sptr_layer = std::dynamic_pointer_cast(r.sptr_shape); - sptr_layer->psi_geom*=(180./sim_pi); - sptr_layer->theta_geom*=(180./sim_pi); - sptr_layer->phi_geom*=(180./sim_pi); + sptr_layer->psi_geom*=(180./simcoon::pi); + sptr_layer->theta_geom*=(180./simcoon::pi); + sptr_layer->phi_geom*=(180./simcoon::pi); paramphases << r.sptr_matprops->number << "\t" << r.sptr_matprops->umat_name << "\t" << r.sptr_matprops->save << "\t" << sptr_layer->concentration << "\t" << r.sptr_matprops->psi_mat << "\t" << r.sptr_matprops->theta_mat << "\t" << r.sptr_matprops->phi_mat << "\t" << sptr_layer->psi_geom << "\t" << sptr_layer->theta_geom << "\t" << sptr_layer->phi_geom << "\t" << r.sptr_matprops->nprops << "\t" << r.sptr_sv_global->nstatev; @@ -105,14 +105,14 @@ void write_ellipsoid(phase_characteristics &rve, const string &path_data, const for(auto r : rve.sub_phases) { - r.sptr_matprops->psi_mat*=(180./sim_pi); - r.sptr_matprops->theta_mat*=(180./sim_pi); - r.sptr_matprops->phi_mat*=(180./sim_pi); + r.sptr_matprops->psi_mat*=(180./simcoon::pi); + r.sptr_matprops->theta_mat*=(180./simcoon::pi); + r.sptr_matprops->phi_mat*=(180./simcoon::pi); auto sptr_ellipsoid = std::dynamic_pointer_cast(r.sptr_shape); - sptr_ellipsoid->psi_geom*=(180./sim_pi); - sptr_ellipsoid->theta_geom*=(180./sim_pi); - sptr_ellipsoid->phi_geom*=(180./sim_pi); + sptr_ellipsoid->psi_geom*=(180./simcoon::pi); + sptr_ellipsoid->theta_geom*=(180./simcoon::pi); + sptr_ellipsoid->phi_geom*=(180./simcoon::pi); paramphases << r.sptr_matprops->number << "\t" << sptr_ellipsoid->coatingof << "\t" << r.sptr_matprops->umat_name << "\t" << r.sptr_matprops->save << "\t" << sptr_ellipsoid->concentration << "\t" << r.sptr_matprops->psi_mat << "\t" << r.sptr_matprops->theta_mat << "\t" << r.sptr_matprops->phi_mat << "\t" << sptr_ellipsoid->a1 << "\t" << sptr_ellipsoid->a2 << "\t" << sptr_ellipsoid->a3 << "\t" << sptr_ellipsoid->psi_geom << "\t" << sptr_ellipsoid->theta_geom << "\t" << sptr_ellipsoid->phi_geom << "\t" << r.sptr_matprops->nprops << "\t" << r.sptr_sv_global->nstatev; @@ -136,14 +136,14 @@ void write_cylinder(phase_characteristics &rve, const string &path_data, const s for(auto r : rve.sub_phases) { - r.sptr_matprops->psi_mat*=(180./sim_pi); - r.sptr_matprops->theta_mat*=(180./sim_pi); - r.sptr_matprops->phi_mat*=(180./sim_pi); + r.sptr_matprops->psi_mat*=(180./simcoon::pi); + r.sptr_matprops->theta_mat*=(180./simcoon::pi); + r.sptr_matprops->phi_mat*=(180./simcoon::pi); auto sptr_cylinder = std::dynamic_pointer_cast(r.sptr_shape); - sptr_cylinder->psi_geom*=(180./sim_pi); - sptr_cylinder->theta_geom*=(180./sim_pi); - sptr_cylinder->phi_geom*=(180./sim_pi); + sptr_cylinder->psi_geom*=(180./simcoon::pi); + sptr_cylinder->theta_geom*=(180./simcoon::pi); + sptr_cylinder->phi_geom*=(180./simcoon::pi); paramphases << r.sptr_matprops->number << "\t" << sptr_cylinder->coatingof << "\t" << r.sptr_matprops->umat_name << "\t" << r.sptr_matprops->save << "\t" << sptr_cylinder->concentration << "\t" << r.sptr_matprops->psi_mat << "\t" << r.sptr_matprops->theta_mat << "\t" << r.sptr_matprops->phi_mat << "\t" << sptr_cylinder->L << "\t" << sptr_cylinder->R << "\t" << sptr_cylinder->psi_geom << "\t" << sptr_cylinder->theta_geom << "\t" << sptr_cylinder->phi_geom << "\t" << r.sptr_matprops->nprops << "\t" << r.sptr_sv_global->nstatev; diff --git a/src/Simulation/Solver/read.cpp b/src/Simulation/Solver/read.cpp index 1c51264bf..da899d44a 100755 --- a/src/Simulation/Solver/read.cpp +++ b/src/Simulation/Solver/read.cpp @@ -132,19 +132,21 @@ void read_matprops(string &umat_name, unsigned int &nprops, vec &props, unsigned props = zeros(nprops); - propsmat.open(path_materialfile, ios::in); - if(!propsmat) { - throw runtime_error("Cannot reopen material file " + materialfile + " in folder " + path_data); - } - + propsmat.open(path_materialfile, ios::in); + if(!propsmat) { + throw runtime_error("Cannot reopen material file " + materialfile + " in folder " + path_data); + } + + // Read orientation angles and properties from material file propsmat >> buffer >> buffer >> buffer >> buffer >> buffer >> buffer >> buffer >> buffer >> buffer >> psi_rve >> buffer >> theta_rve >> buffer >> phi_rve >> buffer; - - for(unsigned int i=0;i> buffer >> props(i); - - psi_rve*=(sim_pi/180.); - theta_rve*=(sim_pi/180.); - phi_rve*=(sim_pi/180.); + } + + psi_rve *= (simcoon::pi / 180.); + theta_rve *= (simcoon::pi / 180.); + phi_rve *= (simcoon::pi / 180.); propsmat.close(); } @@ -272,7 +274,7 @@ void check_path_output(const std::vector &blocks, const solver_output &so } } else if(so.o_type(i) == 2) { - if((fmod(1, so.o_tfreq(i)) > sim_limit)||(fmod(so.o_tfreq(i), sptr_meca->Dn_inc) > 0.)) { + if((fmod(1, so.o_tfreq(i)) > simcoon::limit)||(fmod(so.o_tfreq(i), sptr_meca->Dn_inc) > 0.)) { cout << "The output tfreq is not compatible with the time of increments of the step)" << endl; break; } @@ -303,7 +305,7 @@ void check_path_output(const std::vector &blocks, const solver_output &so } } else if(so.o_type(i) == 2) { - if((fmod(1, so.o_tfreq(i)) > sim_limit)||(fmod(so.o_tfreq(i), sptr_thermomeca->Dn_inc) > 0.)) { + if((fmod(1, so.o_tfreq(i)) > simcoon::limit)||(fmod(so.o_tfreq(i), sptr_thermomeca->Dn_inc) > 0.)) { cout << "The output tfreq is not compatible with the time of increments of the step)" << endl; break; } diff --git a/src/Simulation/Solver/solver.cpp b/src/Simulation/Solver/solver.cpp index a0d3396ea..6f4f075b3 100755 --- a/src/Simulation/Solver/solver.cpp +++ b/src/Simulation/Solver/solver.cpp @@ -306,7 +306,7 @@ void solver(const string &umat_name, const vec &props, const unsigned int &nstat sv_M->DEtot = t2v_strain(Green_Lagrange(sv_M->F1)) - sv_M->Etot; - if (DTime > sim_iota) + if (DTime > simcoon::iota) D = sv_M->Detot/DTime; else D = zeros(3,3); @@ -576,7 +576,7 @@ void solver(const string &umat_name, const vec &props, const unsigned int &nstat if(corate_type == 2) { logarithmic(sv_M->DR, D, Omega, DTime, sv_M->F0, sv_M->F1); } - if (DTime > sim_iota) + if (DTime > simcoon::iota) D = sv_M->Detot/DTime; else D = zeros(3,3); @@ -666,7 +666,7 @@ void solver(const string &umat_name, const vec &props, const unsigned int &nstat error = norm(residual, 2.); if(tnew_dt < 1.) { - if((fabs(Dtinc_cur - sptr_meca->Dn_mini) > sim_iota)||(inforce_solver == 0)) { + if((fabs(Dtinc_cur - sptr_meca->Dn_mini) > simcoon::iota)||(inforce_solver == 0)) { compteur = maxiter_solver; } } @@ -674,7 +674,7 @@ void solver(const string &umat_name, const vec &props, const unsigned int &nstat } } -/* if((fabs(Dtinc_cur - sptr_meca->Dn_mini) < sim_iota)&&(tnew_dt < 1.)) { +/* if((fabs(Dtinc_cur - sptr_meca->Dn_mini) < simcoon::iota)&&(tnew_dt < 1.)) { // cout << "The subroutine has required a step reduction lower than the minimal indicated at" << sptr_meca->number << " inc: " << inc << " and fraction:" << tinc << "\n"; //The solver has been inforced! return; @@ -1003,7 +1003,7 @@ void solver(const string &umat_name, const vec &props, const unsigned int &nstat error = norm(residual, 2.); if(tnew_dt < 1.) { - if((fabs(Dtinc_cur - sptr_thermomeca->Dn_mini) > sim_iota)||(inforce_solver == 0)) { + if((fabs(Dtinc_cur - sptr_thermomeca->Dn_mini) > simcoon::iota)||(inforce_solver == 0)) { compteur = maxiter_solver; } } @@ -1012,7 +1012,7 @@ void solver(const string &umat_name, const vec &props, const unsigned int &nstat } -/* if((fabs(Dtinc_cur - sptr_thermomeca->Dn_mini) < sim_iota)&&(tnew_dt < 1.)) { +/* if((fabs(Dtinc_cur - sptr_thermomeca->Dn_mini) < simcoon::iota)&&(tnew_dt < 1.)) { cout << "The subroutine has required a step reduction lower than the minimal indicated at" << sptr_thermomeca->number << " inc: " << inc << " and fraction:" << tinc << "\n"; //The solver has been inforced! return; diff --git a/src/Simulation/Solver/step_meca.cpp b/src/Simulation/Solver/step_meca.cpp index dbf5ac202..b30b09606 100755 --- a/src/Simulation/Solver/step_meca.cpp +++ b/src/Simulation/Solver/step_meca.cpp @@ -158,7 +158,7 @@ void step_meca::generate(const double &mTime, const vec &mEtot, const vec &msigm if (mode == 2) { double sum_ = 0.; for(int k = 0 ; k < ninc ; k++){ - inc_coef(k) = cos(sim_pi+ (k+1)*2.*sim_pi/(ninc+1))+1.; + inc_coef(k) = cos(simcoon::pi+ (k+1)*2.*simcoon::pi/(ninc+1))+1.; sum_ += inc_coef(k); } inc_coef = inc_coef*ninc/sum_; @@ -307,7 +307,7 @@ void step_meca::generate_kin(const double &mTime, const mat &mF, const double &m if (mode == 2) { double sum_ = 0.; for(int k = 0 ; k < ninc ; k++){ - inc_coef(k) = cos(sim_pi+ (k+1)*2.*sim_pi/(ninc+1))+1.; + inc_coef(k) = cos(simcoon::pi+ (k+1)*2.*simcoon::pi/(ninc+1))+1.; sum_ += inc_coef(k); } inc_coef = inc_coef*ninc/sum_; diff --git a/src/Simulation/Solver/step_thermomeca.cpp b/src/Simulation/Solver/step_thermomeca.cpp index 48f3b1a4d..1f22162c3 100755 --- a/src/Simulation/Solver/step_thermomeca.cpp +++ b/src/Simulation/Solver/step_thermomeca.cpp @@ -155,7 +155,7 @@ void step_thermomeca::generate(const double &mTime, const vec &mEtot, const vec if (mode == 2) { double sum_ = 0.; for(int k = 0 ; k < ninc ; k++){ - inc_coef(k) = cos(sim_pi+ (k+1)*2.*sim_pi/(ninc+1))+1.; + inc_coef(k) = cos(simcoon::pi+ (k+1)*2.*simcoon::pi/(ninc+1))+1.; sum_ += inc_coef(k); } inc_coef = inc_coef*ninc/sum_; @@ -331,7 +331,7 @@ void step_thermomeca::generate_kin(const double &mTime, const mat &mF, const dou if (mode == 2) { double sum_ = 0.; for(int k = 0 ; k < ninc ; k++){ - inc_coef(k) = cos(sim_pi+ (k+1)*2.*sim_pi/(ninc+1))+1.; + inc_coef(k) = cos(simcoon::pi+ (k+1)*2.*simcoon::pi/(ninc+1))+1.; sum_ += inc_coef(k); } inc_coef = inc_coef*ninc/sum_; diff --git a/test/Identification/Tidentification.cpp b/test/Identification/Tidentification.cpp index 6157c1a2a..dd4679439 100755 --- a/test/Identification/Tidentification.cpp +++ b/test/Identification/Tidentification.cpp @@ -101,7 +101,7 @@ TEST(Tidentification, identification_layers) read_layer(rve_layer_0, path_results, "Nlayers0.dat"); for(unsigned int i=0; iprops(i)) > sim_iota) { + if (fabs(rve_layer_0.sptr_matprops->props(i)) > simcoon::iota) { EXPECT_LT( pow(pow(rve_layer_0.sptr_matprops->props(i),2.) - pow(rve_layer_1.sptr_matprops->props(i),2.),0.5)/fabs(rve_layer_0.sptr_matprops->props(i)),1.E-6); } else { diff --git a/test/Libraries/Continuum_mechanics/Tcontimech.cpp b/test/Libraries/Continuum_mechanics/Tcontimech.cpp index 0a555f807..156297c6e 100755 --- a/test/Libraries/Continuum_mechanics/Tcontimech.cpp +++ b/test/Libraries/Continuum_mechanics/Tcontimech.cpp @@ -42,8 +42,8 @@ TEST(Tcontimech, dev_sph) mat testsph = (1./3.)*eye(3,3); - EXPECT_LT(norm(sph(test) - testsph,2),sim_iota); - EXPECT_LT(norm(dev(test) - testdev,2),sim_iota); + EXPECT_LT(norm(sph(test) - testsph,2),simcoon::iota); + EXPECT_LT(norm(dev(test) - testdev,2),simcoon::iota); } TEST(Tcontimech, tr_dev_Mises_eta) @@ -60,8 +60,8 @@ TEST(Tcontimech, tr_dev_Mises_eta) //Test of tr function double tr1 = tr(zeros(6)); double tr2 = tr(test); - EXPECT_LT(fabs(tr1 - 0.),sim_iota); - EXPECT_LT(fabs(tr2 - 12.),sim_iota); + EXPECT_LT(fabs(tr1 - 0.),simcoon::iota); + EXPECT_LT(fabs(tr2 - 12.),simcoon::iota); //Test of dev function vec vide = zeros(6); @@ -75,12 +75,12 @@ TEST(Tcontimech, tr_dev_Mises_eta) testdev(3) = 8.; testdev(4) = 3.; testdev(5) = 7.; - EXPECT_LT(fabs(trdev2 - 0.),sim_iota); - EXPECT_LT(norm(dev2 - testdev,2),sim_iota); + EXPECT_LT(fabs(trdev2 - 0.),simcoon::iota); + EXPECT_LT(norm(dev2 - testdev,2),simcoon::iota); //Test of Mises_stress function double dstress = 3.*sqrt(42.); - EXPECT_LT(fabs(Mises_stress(test) - dstress),sim_iota); + EXPECT_LT(fabs(Mises_stress(test) - dstress),simcoon::iota); //Test of eta_stress function dev1 = eta_stress(zeros(6)); @@ -92,18 +92,18 @@ TEST(Tcontimech, tr_dev_Mises_eta) testeta(3) = 8.*2; testeta(4) = 3.*2; testeta(5) = 7.*2; - EXPECT_LT(norm(dev1,2),sim_iota); - EXPECT_LT(norm(dev2 - (3./2.)*testeta/dstress,2),sim_iota); + EXPECT_LT(norm(dev1,2),simcoon::iota); + EXPECT_LT(norm(dev2 - (3./2.)*testeta/dstress,2),simcoon::iota); //Test of Mises_strain function double dstrain = sqrt(46.); - EXPECT_LT(fabs(Mises_strain(test) - dstrain),sim_iota); + EXPECT_LT(fabs(Mises_strain(test) - dstrain),simcoon::iota); //Test of eta_strain function dev1 = eta_strain(zeros(6)); dev2 = eta_strain(test); - EXPECT_LT(norm(dev1,2),sim_iota); - EXPECT_LT(norm(dev2 - (2./3.)*testdev/dstrain,2),sim_iota); + EXPECT_LT(norm(dev1,2),simcoon::iota); + EXPECT_LT(norm(dev2 - (2./3.)*testdev/dstrain,2),simcoon::iota); } @@ -121,26 +121,26 @@ TEST(Tcontimech, J2_J3) //Test of J2_stress function double J2_stress1 = J2_stress(zeros(6)); double J2_stress2 = J2_stress(test); - EXPECT_LT(fabs(J2_stress1),sim_iota); - EXPECT_LT(fabs(J2_stress2 - 126.),sim_iota); + EXPECT_LT(fabs(J2_stress1),simcoon::iota); + EXPECT_LT(fabs(J2_stress2 - 126.),simcoon::iota); //Test of J2_strain function double J2_strain1 = J2_strain(zeros(6)); double J2_strain2 = J2_strain(test); - EXPECT_LT(fabs(J2_strain1),sim_iota); - EXPECT_LT(fabs(J2_strain2 - 34.5),sim_iota); + EXPECT_LT(fabs(J2_strain1),simcoon::iota); + EXPECT_LT(fabs(J2_strain2 - 34.5),simcoon::iota); //Test of J3_stress function double J3_stress1 = J3_stress(zeros(6)); double J3_stress2 = J3_stress(test); - EXPECT_LT(fabs(J3_stress1),sim_iota); - EXPECT_LT(fabs(J3_stress2 - 226.),sim_iota); + EXPECT_LT(fabs(J3_stress1),simcoon::iota); + EXPECT_LT(fabs(J3_stress2 - 226.),simcoon::iota); //Test of J3_stress function double J3_strain1 = J3_strain(zeros(6)); double J3_strain2 = J3_strain(test); - EXPECT_LT(fabs(J3_strain1),sim_iota); - EXPECT_LT(fabs(J3_strain2 - 14.5),sim_iota); + EXPECT_LT(fabs(J3_strain1),simcoon::iota); + EXPECT_LT(fabs(J3_strain2 - 14.5),simcoon::iota); } @@ -171,21 +171,21 @@ TEST(Tcontimech, ellipsoid) vec sigma_in = {4., 5., 6., 2., 1., 1.5}; vec test_sig_int = {5.,2.5}; vec normal = normal_ellipsoid(u,v,a1,a2,a3); - EXPECT_LT(norm(normal-test,2),sim_iota); + EXPECT_LT(norm(normal-test,2),simcoon::iota); u = 0.; - v =sim_pi/2.; + v =simcoon::pi/2.; test = { 1, 0, 0 }; normal = normal_ellipsoid(u,v,a1,a2,a3); - EXPECT_LT(norm(normal-test,2),sim_iota); + EXPECT_LT(norm(normal-test,2),simcoon::iota); - u =sim_pi/2.; - v =sim_pi/2.; + u =simcoon::pi/2.; + v =simcoon::pi/2.; test = { 0, 1, 0 }; normal = normal_ellipsoid(u,v,a1,a2,a3); - EXPECT_LT(norm(normal-test,2),sim_iota); + EXPECT_LT(norm(normal-test,2),simcoon::iota); vec sig_int = sigma_int(sigma_in,u,v,a1,a2,a3); - EXPECT_LT(norm(sig_int-test_sig_int,2),sim_iota); + EXPECT_LT(norm(sig_int-test_sig_int,2),simcoon::iota); } TEST(Tcontimech, P_ijkl) @@ -260,7 +260,7 @@ TEST(Tcontimech, P_ijkl) mat pikjla = p_ikjl(a); mat pikjlb = p_ikjl(b); - EXPECT_LT(norm(pikjla - Ireal,2),sim_iota); - EXPECT_LT(norm(pikjlb - ones(6,6),2),sim_iota); - EXPECT_LT(norm(pikjlc - result,2),sim_iota); + EXPECT_LT(norm(pikjla - Ireal,2),simcoon::iota); + EXPECT_LT(norm(pikjlb - ones(6,6),2),simcoon::iota); + EXPECT_LT(norm(pikjlc - result,2),simcoon::iota); } diff --git a/test/Libraries/Continuum_mechanics/Tcriteria.cpp b/test/Libraries/Continuum_mechanics/Tcriteria.cpp index cd753df25..c1fef5992 100755 --- a/test/Libraries/Continuum_mechanics/Tcriteria.cpp +++ b/test/Libraries/Continuum_mechanics/Tcriteria.cpp @@ -38,9 +38,9 @@ TEST(Tcriteria, Aniso) // Check that Tresca is equal to 400 in that case double test_tresca = Tresca_stress(sigma); - EXPECT_LT(test_tresca - 400., sim_iota); + EXPECT_LT(test_tresca - 400., simcoon::iota); // Check that Drucker is equal to 400 in that case double test_prager = Drucker_stress(sigma, b, n); - EXPECT_LT(test_prager - 400., sim_iota); + EXPECT_LT(test_prager - 400., simcoon::iota); } \ No newline at end of file diff --git a/test/Libraries/Continuum_mechanics/Thyperelastic.cpp b/test/Libraries/Continuum_mechanics/Thyperelastic.cpp index 6616a0da8..d7622c618 100755 --- a/test/Libraries/Continuum_mechanics/Thyperelastic.cpp +++ b/test/Libraries/Continuum_mechanics/Thyperelastic.cpp @@ -21,8 +21,8 @@ #include #include - #include +#include #include #include @@ -51,12 +51,25 @@ TEST(Thyperelastic, isochoric_invariants) mat b_rand = simcoon::v2t_strain(randu(6))+eye(3,3); - double J = sqrt(det(b_rand)); + double J; + try { + J = sqrt(det(b_rand)); + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside Thyperelastic."); + } mat b_bar = pow(J,-2./3.)*b_rand; vec I = zeros(3); + I(0) = trace(b_bar); - I(1) = 0.5*(pow(trace(b_bar),2.)-trace(powmat(b_bar,2))); + mat b_bar2; + try { + I(1) = 0.5*(pow(trace(b_bar),2.)-trace(powmat(b_bar,2))); + } catch (const std::runtime_error &e) { + cerr << "Error in det: " << e.what() << endl; + throw simcoon::exception_powmat("Error in powmat function inside Thyperelastic."); + } I(2) = 1.; vec I_rand = simcoon::isochoric_invariants(b_rand); @@ -96,7 +109,14 @@ TEST(Thyperelastic, isochoric_pstretch) mat V_rand = simcoon::v2t_strain(randu(6)) + eye(3,3); mat b_rand = powmat(V_rand,2); - double J = det(V_rand); + + double J; + try { + J = det(V_rand); + } catch (const std::runtime_error &e) { + cerr << "Error in det or inv: throw inv exception " << e.what() << endl; + throw simcoon::exception_det("Error in det function inside Thyperelastic."); + } vec lambda = eig_sym(V_rand); vec lambda_bar = pow(J,-1./3.)*lambda; diff --git a/test/Libraries/Continuum_mechanics/Tobjective_rates.cpp b/test/Libraries/Continuum_mechanics/Tobjective_rates.cpp index 5bf24ed37..09ae375b3 100755 --- a/test/Libraries/Continuum_mechanics/Tobjective_rates.cpp +++ b/test/Libraries/Continuum_mechanics/Tobjective_rates.cpp @@ -67,7 +67,7 @@ TEST(Tobjective_rates, get_B) double f_z = 0.; for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { - if ((i!=j)&&(fabs(bi(i)-bi(j))>sim_iota)) { + if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) { Bij = Bi.col(i)*(Bi.col(j)).t(); Bij_ = mat_FTensor2(Bij); f_z = (1.+(bi(i)/bi(j)))/(1.-(bi(i)/bi(j)))+2./log(bi(i)/bi(j)); @@ -140,7 +140,7 @@ TEST(Tobjective_rates, logarithmic_functions) BBBB_(k,l,m,n) = Bij_(k,l)*Bij_(m,n); for (unsigned int i=0; i<3; i++) { for (unsigned int j=0; j<3; j++) { - if ((i!=j)&&(fabs(bi(i)-bi(j))>sim_iota)) { + if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) { Bij = Bi.col(i)*(Bi.col(j)).t(); Bij_ = mat_FTensor2(Bij); f_z = (1.+(bi(i)/bi(j)))/(1.-(bi(i)/bi(j)))+2./log(bi(i)/bi(j)); diff --git a/test/Libraries/Continuum_mechanics/Ttransfer.cpp b/test/Libraries/Continuum_mechanics/Ttransfer.cpp index b424c7ba2..9411fc27c 100755 --- a/test/Libraries/Continuum_mechanics/Ttransfer.cpp +++ b/test/Libraries/Continuum_mechanics/Ttransfer.cpp @@ -69,27 +69,27 @@ TEST(Ttransfer, v2t_func) //Test of v2t_strain function strain = v2t_strain(zeros(6)); - EXPECT_LT(norm(strain,2),sim_iota); + EXPECT_LT(norm(strain,2),simcoon::iota); strain = v2t_strain(test); - EXPECT_LT(norm(strain - teststrain,2),sim_iota); + EXPECT_LT(norm(strain - teststrain,2),simcoon::iota); //Test of t2v_strain function vec dev1 = t2v_strain(zeros(3,3)); vec dev2 = t2v_strain(teststrain); - EXPECT_LT(norm(dev1,2),sim_iota); - EXPECT_LT(norm(dev2 - test,2),sim_iota); + EXPECT_LT(norm(dev1,2),simcoon::iota); + EXPECT_LT(norm(dev2 - test,2),simcoon::iota); //Test of v2t_stress function stress = v2t_stress(zeros(6)); - EXPECT_LT(norm(stress,2),sim_iota); + EXPECT_LT(norm(stress,2),simcoon::iota); stress = v2t_stress(test); - EXPECT_LT(norm(stress - teststress,2),sim_iota); + EXPECT_LT(norm(stress - teststress,2),simcoon::iota); //Test of t2v_stress function dev1 = t2v_stress(zeros(3,3)); dev2 = t2v_stress(teststress); - EXPECT_LT(norm(dev1,2),sim_iota); - EXPECT_LT(norm(dev2 - test,2),sim_iota); + EXPECT_LT(norm(dev1,2),simcoon::iota); + EXPECT_LT(norm(dev2 - test,2),simcoon::iota); } @@ -174,9 +174,9 @@ TEST(Ttransfer, FTensor_transfer) C = mat_FTensor4(L); mat L2 = FTensor4_mat(C); - EXPECT_LT(norm(testmat2 - testmat,2),sim_iota); - EXPECT_LT(norm(test_strain - test,2),sim_iota); - EXPECT_LT(norm(test_stress - test,2),sim_iota); - EXPECT_LT(norm(L2 - L,2),sim_iota); + EXPECT_LT(norm(testmat2 - testmat,2),simcoon::iota); + EXPECT_LT(norm(test_strain - test,2),simcoon::iota); + EXPECT_LT(norm(test_stress - test,2),simcoon::iota); + EXPECT_LT(norm(L2 - L,2),simcoon::iota); } diff --git a/test/Libraries/Maths/Trotation.cpp b/test/Libraries/Maths/Trotation.cpp index 4659fc780..fb108be19 100755 --- a/test/Libraries/Maths/Trotation.cpp +++ b/test/Libraries/Maths/Trotation.cpp @@ -38,9 +38,9 @@ TEST(Trotation, basic_rotation) double theta = 42.; double phi = 165; - psi *= (sim_pi/180.); - theta *= (sim_pi/180.); - phi *= (sim_pi/180.); + psi *= (simcoon::pi/180.); + theta *= (simcoon::pi/180.); + phi *= (simcoon::pi/180.); mat R1 = { {cos(psi),-sin(psi),0}, {sin(psi),cos(psi), 0}, {0,0,1}}; mat R2 = { {1,0,0}, {0,cos(theta),-sin(theta)}, {0,sin(theta),cos(theta)} }; @@ -93,9 +93,9 @@ TEST(Trotation, rotation) double theta = 32.; double phi = -4.5; - psi *= (sim_pi/180.); - theta *= (sim_pi/180.); - phi *= (sim_pi/180.); + psi *= (simcoon::pi/180.); + theta *= (simcoon::pi/180.); + phi *= (simcoon::pi/180.); //Rotation mat R1 = { {cos(psi),-sin(psi),0}, {sin(psi),cos(psi), 0}, {0,0,1}}; diff --git a/test/Umats/LOG_int/TLOG_int.cpp b/test/Umats/LOG_int/TLOG_int.cpp index 8889b63a2..42499a96f 100755 --- a/test/Umats/LOG_int/TLOG_int.cpp +++ b/test/Umats/LOG_int/TLOG_int.cpp @@ -28,6 +28,7 @@ #include #include +#include #include #include #include @@ -94,8 +95,15 @@ TEST(TLOG_int, TLOG_int_solver) F_test(2,0) = 0.; F_test(2,1) = 0.; F_test(2,2) = 1.; - - vec e_tot_log_test = t2v_strain(0.5*logmat_sympd(L_Cauchy_Green(F_test))); + + vec e_tot_log_test; + try { + e_tot_log_test = t2v_strain(0.5*logmat_sympd(L_Cauchy_Green(F_test))); + } catch (const std::runtime_error &e) { + cerr << "Error in logmat_sympd: " << e.what() << endl; + throw simcoon::exception_logmat_sympd("Error in logmat_sympd function inside TLOG_int."); + } + unsigned int n_rows_results = R.n_rows; vec e_tot_log = zeros(6); diff --git a/testBin/Umats/UMABA/external/UMAT_ABAQUS_ELASTIC.for b/testBin/Umats/UMABA/external/UMAT_ABAQUS_ELASTIC.f similarity index 100% rename from testBin/Umats/UMABA/external/UMAT_ABAQUS_ELASTIC.for rename to testBin/Umats/UMABA/external/UMAT_ABAQUS_ELASTIC.f