From 8cf72e2233b3208b287c014988f36a84ca0f291c Mon Sep 17 00:00:00 2001 From: Thomas Helfer Date: Wed, 6 Sep 2023 10:50:41 +0200 Subject: [PATCH] First implementation. Update of the bindings --- .../MGIS/Behaviour/MaterialStateManager.h | 38 +++++++ bindings/c/include/MGIS/Behaviour/State.h | 7 ++ bindings/c/src/MaterialStateManager.cxx | 72 +++++++++++- bindings/c/src/State.cxx | 9 ++ bindings/fortran/src/mgis_behaviour.f95 | 103 ++++++++++++++++++ bindings/python/src/MaterialStateManager.cxx | 19 ++++ bindings/python/src/State.cxx | 1 + docs/web/CMakeLists.txt | 1 + docs/web/mgis-template.html | 3 +- docs/web/release-notes-2.2.md | 64 +++++++++++ .../MGIS/Behaviour/MaterialDataManager.hxx | 12 +- .../MGIS/Behaviour/MaterialStateManager.hxx | 41 ++++++- include/MGIS/Behaviour/State.hxx | 2 + src/Integrate.cxx | 67 ++++++++++-- src/MaterialStateManager.cxx | 37 +++++++ src/State.cxx | 1 + tests/CMakeLists.txt | 37 ++++++- tests/ComputeSpeedOfSound.mfront | 21 ++++ tests/ComputeSpeedOfSoundTest.cxx | 63 +++++++++++ tests/ComputeSpeedOfSoundTest2.cxx | 62 +++++++++++ 20 files changed, 639 insertions(+), 21 deletions(-) create mode 100644 docs/web/release-notes-2.2.md create mode 100644 tests/ComputeSpeedOfSound.mfront create mode 100644 tests/ComputeSpeedOfSoundTest.cxx create mode 100644 tests/ComputeSpeedOfSoundTest2.cxx diff --git a/bindings/c/include/MGIS/Behaviour/MaterialStateManager.h b/bindings/c/include/MGIS/Behaviour/MaterialStateManager.h index 45ca20790..99e8943b4 100644 --- a/bindings/c/include/MGIS/Behaviour/MaterialStateManager.h +++ b/bindings/c/include/MGIS/Behaviour/MaterialStateManager.h @@ -198,6 +198,44 @@ MGIS_C_EXPORT mgis_status mgis_bv_material_state_manager_get_non_uniform_material_property( mgis_real** const, mgis_bv_MaterialStateManager* const, const char* const); +/*! + * \brief set the value of an uniform scalar mass density + * \param[in] s: state manager + * \param[in] v: value of the mass density + */ +MGIS_C_EXPORT mgis_status +mgis_bv_material_state_manager_set_uniform_scalar_mass_density( + mgis_bv_MaterialStateManager* const, const mgis_real); +/*! + * \brief set the value of an uniform mass density + * \param[in] s: state manager + * \param[in] v: values of the mass density + * \param[in] sm: storage mode + */ +MGIS_C_EXPORT mgis_status +mgis_bv_material_state_manager_set_non_uniform_mass_density( + mgis_bv_MaterialStateManager* const, + mgis_real* const, + const mgis_bv_MaterialStateManagerStorageMode); + +MGIS_C_EXPORT mgis_status +mgis_bv_material_state_manager_is_mass_density_defined( + int* const, const mgis_bv_MaterialStateManager* const); + +MGIS_C_EXPORT mgis_status +mgis_bv_material_state_manager_is_mass_density_uniform( + int* const, const mgis_bv_MaterialStateManager* const); + +MGIS_C_EXPORT mgis_status +mgis_bv_material_state_manager_get_uniform_mass_density( + mgis_real* const, mgis_bv_MaterialStateManager* const); + +MGIS_C_EXPORT mgis_status +mgis_bv_material_state_manager_get_non_uniform_mass_density( + mgis_real** const, mgis_bv_MaterialStateManager* const); + +// + MGIS_C_EXPORT mgis_status mgis_bv_material_state_manager_set_uniform_scalar_external_state_variable( mgis_bv_MaterialStateManager* const, const char* const, const mgis_real); diff --git a/bindings/c/include/MGIS/Behaviour/State.h b/bindings/c/include/MGIS/Behaviour/State.h index 0d2344416..96c477069 100644 --- a/bindings/c/include/MGIS/Behaviour/State.h +++ b/bindings/c/include/MGIS/Behaviour/State.h @@ -34,6 +34,13 @@ using mgis_bv_State = mgis::behaviour::State; typedef struct mgis_bv_State mgis_bv_State; #endif /* __cplusplus */ +/*! + * \brief set the mass density + * \param[out] s: state + * \param[in] v: value + */ +MGIS_C_EXPORT mgis_status mgis_bv_state_set_mass_density(mgis_bv_State* const, + const mgis_real); /*! * \brief set a gradient' value in a state * \param[out] s: state diff --git a/bindings/c/src/MaterialStateManager.cxx b/bindings/c/src/MaterialStateManager.cxx index 289c0bdb0..9dacd0268 100644 --- a/bindings/c/src/MaterialStateManager.cxx +++ b/bindings/c/src/MaterialStateManager.cxx @@ -266,7 +266,7 @@ mgis_status mgis_bv_material_state_manager_set_non_uniform_material_property( return mgis_report_failure("null state manager"); } if (v == nullptr) { - return mgis_report_failure("null state values"); + return mgis_report_failure("null material property values"); } try { const auto& mp = getVariable(m->b.mps, n); @@ -317,6 +317,76 @@ mgis_status mgis_bv_material_state_manager_is_material_property_uniform( return mgis_report_success(); } // end of mgis_bv_material_state_manager_is_material_property_uniform +MGIS_C_EXPORT mgis_status +mgis_bv_material_state_manager_set_uniform_scalar_mass_density( + mgis_bv_MaterialStateManager* const m, + const mgis_real v) { + if (m == nullptr) { + return mgis_report_failure("null state manager"); + } + try { + setMassDensity(*m, v); + } catch (...) { + return mgis_handle_cxx_exception(); + } + return mgis_report_success(); +} // end of mgis_bv_material_state_manager_set_uniform_scalar_mass_density + +mgis_status mgis_bv_material_state_manager_set_non_uniform_mass_density( + mgis_bv_MaterialStateManager* const m, + mgis_real* const v, + const mgis_bv_MaterialStateManagerStorageMode s) { + using index_type = mgis::span::index_type; + if (m == nullptr) { + return mgis_report_failure("null state manager"); + } + if (v == nullptr) { + return mgis_report_failure("null state values"); + } + try { + if (s == MGIS_BV_LOCAL_STORAGE) { + setMassDensity(*m, {v, static_cast(m->n)}, + mgis::behaviour::MaterialStateManager::LOCAL_STORAGE); + } else { + setMassDensity(*m, {v, static_cast(m->n)}, + mgis::behaviour::MaterialStateManager::EXTERNAL_STORAGE); + } + } catch (...) { + return mgis_handle_cxx_exception(); + } + return mgis_report_success(); +} // end of mgis_bv_material_state_manager_set_non_uniform_mass_density + +mgis_status mgis_bv_material_state_manager_is_mass_density_defined( + int* const b, + const mgis_bv_MaterialStateManager* const m) { + *b = 0; + if (m == nullptr) { + return mgis_report_failure("null state manager"); + } + try { + *b = isMassDensityDefined(*m); + } catch (...) { + return mgis_handle_cxx_exception(); + } + return mgis_report_success(); +} // end of mgis_bv_material_state_manager_is_mass_density_defined + +mgis_status mgis_bv_material_state_manager_is_mass_density_uniform( + int* const b, + const mgis_bv_MaterialStateManager* const m) { + *b = 0; + if (m == nullptr) { + return mgis_report_failure("null state manager"); + } + try { + *b = isMassDensityUniform(*m); + } catch (...) { + return mgis_handle_cxx_exception(); + } + return mgis_report_success(); +} // end of mgis_bv_material_state_manager_is_mass_density_uniform + MGIS_C_EXPORT mgis_status mgis_bv_material_state_manager_set_uniform_scalar_external_state_variable( mgis_bv_MaterialStateManager* const m, diff --git a/bindings/c/src/State.cxx b/bindings/c/src/State.cxx index eafda7ac2..a449c3c06 100644 --- a/bindings/c/src/State.cxx +++ b/bindings/c/src/State.cxx @@ -18,6 +18,15 @@ extern "C" { +mgis_status mgis_bv_state_set_mass_density(mgis_bv_State* const s, + const mgis_real v) { + if (s == nullptr) { + return mgis_report_failure("invalid argument (null state)"); + } + s->mass_density = v; + return mgis_report_success(); +} // end of mgis_bv_state_set_mass_density + mgis_status mgis_bv_state_set_gradient_by_name(mgis_bv_State* const s, const char* const n, const mgis_real* const v) { diff --git a/bindings/fortran/src/mgis_behaviour.f95 b/bindings/fortran/src/mgis_behaviour.f95 index fb02e88f5..049bb29f4 100644 --- a/bindings/fortran/src/mgis_behaviour.f95 +++ b/bindings/fortran/src/mgis_behaviour.f95 @@ -2111,6 +2111,28 @@ end function free_behaviour_data_wrapper end if end function free_behaviour_data ! + function state_set_mass_density(s, v) result(r) + use mgis_fortran_utilities + use mgis, only: mgis_status + implicit none + interface + function state_set_mass_density_wrapper(s, v) & + bind(c,name = 'mgis_bv_state_set_mass_density') & + result(r) + use, intrinsic :: iso_c_binding, only: c_ptr, c_double + use mgis, only: mgis_status + implicit none + type(c_ptr), intent(in),value :: s + real(c_double), intent(in),value :: v + type(mgis_status) :: r + end function state_set_mass_density_wrapper + end interface + type(State), intent(in) :: s + real(kind=8) :: v + type(mgis_status) :: r + r = state_set_mass_density_wrapper(s%ptr, v) + end function state_set_mass_density + ! function state_set_gradient_by_name(s, n, v) result(r) use, intrinsic :: iso_c_binding, only: c_loc use mgis_fortran_utilities @@ -3130,6 +3152,87 @@ end function msm_is_material_property_uniform end if end function material_state_manager_is_material_property_uniform ! + function material_state_manager_set_uniform_scalar_mass_density(s, v) & + result(r) + use mgis, only: mgis_status + use mgis_fortran_utilities + implicit none + interface + function msm_set_uniform_scalar_mass_density_wrapper(s, v) & + bind(c,name = 'mgis_bv_material_state_manager_set_uniform_scalar_mass_density') & + result(r) + use, intrinsic :: iso_c_binding, only: c_ptr, c_double, c_char + use mgis, only: mgis_status + implicit none + type(c_ptr), intent(in),value :: s + real(kind=c_double), intent(in), value :: v + type(mgis_status) :: r + end function msm_set_uniform_scalar_mass_density_wrapper + end interface + type(MaterialStateManager), intent(in) :: s + real(kind=8), intent(in) :: v + type(mgis_status) :: r + r = msm_set_uniform_scalar_mass_density_wrapper(s %ptr, v) + end function material_state_manager_set_uniform_scalar_mass_density + ! + function material_state_manager_is_mass_density_defined(b, s)& + result(r) + use mgis_fortran_utilities + use mgis, only: mgis_status, MGIS_SUCCESS + implicit none + interface + function msm_is_mass_density_defined(b, s) & + bind(c,name = 'mgis_bv_material_state_manager_is_mass_density_defined') & + result(r) + use, intrinsic :: iso_c_binding, only: c_ptr, c_int, c_char + use mgis, only: mgis_status + implicit none + integer(kind=c_int), intent(out) :: b + type(c_ptr), intent(in),value :: s + type(mgis_status) :: r + end function msm_is_mass_density_defined + end interface + logical, intent(out) :: b + type(MaterialStateManager), intent(in) :: s + type(mgis_status) :: r + integer bc + r = msm_is_mass_density_defined(bc, s%ptr) + if( r % exit_status .eq. MGIS_SUCCESS) then + b = bc .eq. 1 + else + b = .false. + end if + end function material_state_manager_is_mass_density_defined + ! + function material_state_manager_is_mass_density_uniform(b, s)& + result(r) + use mgis_fortran_utilities + use mgis, only: mgis_status, MGIS_SUCCESS + implicit none + interface + function msm_is_mass_density_uniform(b, s) & + bind(c,name = 'mgis_bv_material_state_manager_is_mass_density_uniform') & + result(r) + use, intrinsic :: iso_c_binding, only: c_ptr, c_int, c_char + use mgis, only: mgis_status + implicit none + integer(kind=c_int), intent(out) :: b + type(c_ptr), intent(in),value :: s + type(mgis_status) :: r + end function msm_is_mass_density_uniform + end interface + logical, intent(out) :: b + type(MaterialStateManager), intent(in) :: s + type(mgis_status) :: r + integer bc + r = msm_is_mass_density_uniform(bc, s%ptr) + if( r % exit_status .eq. MGIS_SUCCESS) then + b = bc .eq. 1 + else + b = .false. + end if + end function material_state_manager_is_mass_density_uniform + ! function material_state_manager_get_internal_state_variables_stride(n_isvs, s) & result(r) use, intrinsic :: iso_c_binding, only: c_size_t diff --git a/bindings/python/src/MaterialStateManager.cxx b/bindings/python/src/MaterialStateManager.cxx index a7630c0ec..c0863578f 100644 --- a/bindings/python/src/MaterialStateManager.cxx +++ b/bindings/python/src/MaterialStateManager.cxx @@ -102,6 +102,19 @@ static void MaterialStateManager_setMaterialProperty2( setMaterialProperty(sm, n, mgis::python::mgis_convert_to_span(o), s); } // end of MaterialStateManager_setMaterialProperty +static void MaterialStateManager_setMassDensity( + mgis::behaviour::MaterialStateManager& s, + const mgis::real v) { + mgis::behaviour::setMassDensity(s, v); +} // end of MaterialStateManager_setMassDensity + +static void MaterialStateManager_setMassDensity2( + mgis::behaviour::MaterialStateManager& sm, + const boost::python::object& o, + const mgis::behaviour::MaterialStateManager::StorageMode s) { + setMassDensity(sm, mgis::python::mgis_convert_to_span(o), s); +} // end of MaterialStateManager_setMassDensity + static void MaterialStateManager_setExternalStateVariable( mgis::behaviour::MaterialStateManager& s, const std::string& n, @@ -174,6 +187,8 @@ void declareMaterialStateManager() { &MaterialStateManager_getInternalStateVariables) .def("setMaterialProperty", &MaterialStateManager_setMaterialProperty) .def("setMaterialProperty", &MaterialStateManager_setMaterialProperty2) + .def("setMassDensity", &MaterialStateManager_setMassDensity) + .def("setMassDensity", &MaterialStateManager_setMassDensity2) .def("setExternalStateVariable", &MaterialStateManager_setExternalStateVariable) .def("setExternalStateVariable", @@ -183,6 +198,10 @@ void declareMaterialStateManager() { &MaterialStateManager_setMaterialProperty); boost::python::def("setMaterialProperty", &MaterialStateManager_setMaterialProperty2); + boost::python::def("setMassDensity", + &MaterialStateManager_setMassDensity); + boost::python::def("setMassDensity", + &MaterialStateManager_setMassDensity2); boost::python::def("setExternalStateVariable", &MaterialStateManager_setExternalStateVariable); boost::python::def("setExternalStateVariable", diff --git a/bindings/python/src/State.cxx b/bindings/python/src/State.cxx index 86ee18efe..4cecd1fcb 100644 --- a/bindings/python/src/State.cxx +++ b/bindings/python/src/State.cxx @@ -71,6 +71,7 @@ void declareState() { using mgis::behaviour::Behaviour; using mgis::behaviour::State; boost::python::class_("State", boost::python::no_init) + .add_property("mass_density", &State::mass_density) .add_property("stored_energy", &State::stored_energy) .add_property("dissipated_energy", &State::dissipated_energy) .add_property("gradients", &State_getGradients) diff --git a/docs/web/CMakeLists.txt b/docs/web/CMakeLists.txt index 41132ec3f..6c5014f12 100644 --- a/docs/web/CMakeLists.txt +++ b/docs/web/CMakeLists.txt @@ -71,6 +71,7 @@ if(MGIS_HAVE_PANDOC) mgis_pandoc_generate_html_page(release-notes-1.2.2 "--number-sections" "--toc" "--toc-depth=3") mgis_pandoc_generate_html_page(release-notes-2.0 "--toc" "--toc-depth=3") mgis_pandoc_generate_html_page(release-notes-2.1 "--toc" "--toc-depth=3") + mgis_pandoc_generate_html_page(release-notes-2.2 "--toc" "--toc-depth=3") mgis_pandoc_generate_html_page(orthotropic-behaviours "--toc" "--toc-depth=3") mgis_pandoc_generate_html_page(FEniCSBindings) mgis_pandoc_generate_html_page(mgis_fenics "--toc" "--toc-depth=3") diff --git a/docs/web/mgis-template.html b/docs/web/mgis-template.html index 86430ceff..c9b5c34c6 100644 --- a/docs/web/mgis-template.html +++ b/docs/web/mgis-template.html @@ -89,7 +89,8 @@
  • Version 2.0.x
    • diff --git a/docs/web/release-notes-2.2.md b/docs/web/release-notes-2.2.md new file mode 100644 index 000000000..3ac41a170 --- /dev/null +++ b/docs/web/release-notes-2.2.md @@ -0,0 +1,64 @@ +--- +title: MFrontGenericInterfaceSupport Version 2.2 +author: Thomas Helfer +date: 2020 +lang: en-EN +numbersections: true +documentclass: article +from: markdown+tex_math_single_backslash +geometry: + - margin=2cm +papersize: a4 +link-citations: true +colorlinks: true +figPrefixTemplate: "$$i$$" +tabPrefixTemplate: "$$i$$" +secPrefixTemplate: "$$i$$" +eqnPrefixTemplate: "($$i$$)" +bibliography: bibliography.bib +--- + +# New features + +## Declaration of the mass density in the `MaterialStateManager` class and associated free functions {#sec:mgis:2.2:mass_density} + +The `MaterialStateManager` class now have a data member named `mass_density`. + +The following free functions are available in `C++`: + +- `setMassDensity` which sets the mass density using either a scalar + value (uniform case) or a set of values (non uniform case). +- `isMassDensityDefined` which returns if the mass density has been + defined. +- `isMassDensityUniform` which returns if the mass density is uniform. + This function throws if the mass density has not been defined. + +### `python` bindings + +The `MaterialStateManager` class now have a `setMassDensity` member. A +corresponding free function is also available. + +### `C` bindings + +The following functions are now available: +`mgis_bv_material_state_manager_set_uniform_scalar_mass_density`, +`mgis_bv_material_state_manager_set_non_uniform_mass_density`, +`mgis_bv_material_state_manager_is_mass_density_defined`, +`mgis_bv_material_state_manager_is_mass_density_uniform`, +`mgis_bv_material_state_manager_get_uniform_mass_density` and +`mgis_bv_material_state_manager_get_non_uniform_mass_density`. + +### `fortran` bindings + +The following functions are now available: +`material_state_manager_set_uniform_scalar_mass_density`, +`material_state_manager_is_mass_density_defined` and +`material_state_manager_is_mass_density_uniform`. + +# Issues solved + +## Issue #121: Pass the mass density from the `MaterialDataManager` + +This feature is described in Section @sec:mgis:2.2:mass_density. + +For more details, see . diff --git a/include/MGIS/Behaviour/MaterialDataManager.hxx b/include/MGIS/Behaviour/MaterialDataManager.hxx index c7c6644f8..15b26d280 100644 --- a/include/MGIS/Behaviour/MaterialDataManager.hxx +++ b/include/MGIS/Behaviour/MaterialDataManager.hxx @@ -47,14 +47,18 @@ namespace mgis::behaviour { ~BehaviourIntegrationWorkSpace(); //! \brief a buffer to hold error messages std::vector error_message; - //! material properties at the beginning of the time step + //! \brief material properties at the beginning of the time step std::vector mps0; - //! material properties at the end of the time step + //! \brief material properties at the end of the time step std::vector mps1; - //! external state variables at the beginning of the time step + //! \brief external state variables at the beginning of the time step std::vector esvs0; - //! external state variables at the end of the time step + //! \brief external state variables at the end of the time step std::vector esvs1; + //! \brief mass density at the beginning of the time step + mgis::real rho0; + //! \brief mass density at the end of the time step + mgis::real rho1; }; // end of struct BehaviourIntegrationWorkSpace /*! diff --git a/include/MGIS/Behaviour/MaterialStateManager.hxx b/include/MGIS/Behaviour/MaterialStateManager.hxx index 709cc5a60..b4909ce11 100644 --- a/include/MGIS/Behaviour/MaterialStateManager.hxx +++ b/include/MGIS/Behaviour/MaterialStateManager.hxx @@ -19,6 +19,7 @@ #include #include #include +#include #include "MGIS/Config.hxx" #include "MGIS/Span.hxx" #include "MGIS/StorageMode.hxx" @@ -141,6 +142,14 @@ namespace mgis::behaviour { * case). */ std::map material_properties; + /*! \brief mass density + * The mass density can be uniform or not. + * In the non uniform case, the data can be hold by the structure + * (std::vector) or simply borrow a reference + * (mgis::span + * case). + */ + std::optional mass_density; //! \brief view to the values of the internal state variables mgis::span internal_state_variables; /*! @@ -215,14 +224,38 @@ namespace mgis::behaviour { MGIS_EXPORT bool isMaterialPropertyDefined(const MaterialStateManager&, const mgis::string_view&); /*! - * \brief set the given material property + * \brief chek if the given material property is uniform * \param[out] m: material data manager * \param[in] n: name - * \param[in] v: values - * \param[in] s: storage mode */ MGIS_EXPORT bool isMaterialPropertyUniform(const MaterialStateManager&, const mgis::string_view&); + /*! + * \brief set the mass density + * \param[out] m: material data manager + * \param[in] v: value + */ + MGIS_EXPORT void setMassDensity(MaterialStateManager&, const real); + /*! + * \brief set the mass density + * \param[out] m: material data manager + * \param[in] v: values + * \param[in] s: storage mode + */ + MGIS_EXPORT void setMassDensity(MaterialStateManager&, + const mgis::span&, + const MaterialStateManager::StorageMode = + MaterialStateManager::LOCAL_STORAGE); + /*! + * \return true if the given external state variable is defined. + * \param[out] m: material data manager + */ + MGIS_EXPORT bool isMassDensityDefined(const MaterialStateManager&); + /*! + * \return true if the mass density is uniform + * \param[out] m: material data manager + */ + MGIS_EXPORT bool isMassDensityUniform(const MaterialStateManager&); /*! * \brief set the given external state variable * \param[out] m: material data manager @@ -258,8 +291,6 @@ namespace mgis::behaviour { * \return true if the given external state variable is uniform. * \param[out] m: material data manager * \param[in] n: name - * \param[in] v: values - * \param[in] s: storage mode */ MGIS_EXPORT bool isExternalStateVariableUniform(const MaterialStateManager&, const mgis::string_view&); diff --git a/include/MGIS/Behaviour/State.hxx b/include/MGIS/Behaviour/State.hxx index d3974fe9f..89be1653b 100644 --- a/include/MGIS/Behaviour/State.hxx +++ b/include/MGIS/Behaviour/State.hxx @@ -56,6 +56,8 @@ namespace mgis::behaviour { * files) This output is optional. */ real dissipated_energy; + //! + real mass_density = 0; //! \brief value of the gradients std::vector gradients; //! \brief values of the thermodynamic_forces diff --git a/src/Integrate.cxx b/src/Integrate.cxx index dbee1f953..dd2b999cb 100644 --- a/src/Integrate.cxx +++ b/src/Integrate.cxx @@ -52,13 +52,13 @@ namespace mgis::behaviour::internals { * \brief uniform values are treated immediatly. For spatially variable * fields, we return the information needed to evaluate them */ - static inline std::vector buildEvaluator( + static std::vector buildEvaluators( std::vector& v, std::map& values, const MaterialDataManager& m, const std::vector& ds) { mgis::raise_if(v.size() != getArraySize(ds, m.b.hypothesis), - "buildEvaluator: ill allocated memory"); + "buildEvaluators: ill allocated memory"); // evaluators std::vector evaluators; auto offset = mgis::size_type{}; @@ -66,7 +66,7 @@ namespace mgis::behaviour::internals { const auto s = getVariableSize(d, m.b.hypothesis); auto p = values.find(d.name); if (p == values.end()) { - auto msg = std::string{"buildEvaluator: no variable named '" + d.name + + auto msg = std::string{"buildEvaluators: no variable named '" + d.name + "' declared"}; if (!values.empty()) { msg += "\nThe following variables were declared: "; @@ -91,7 +91,7 @@ namespace mgis::behaviour::internals { if (std::holds_alternative(p->second)) { if (d.type != Variable::SCALAR) { mgis::raise( - "buildEvaluator: invalid type for " + "buildEvaluators: invalid type for " "variable '" + d.name + "'"); } @@ -104,7 +104,36 @@ namespace mgis::behaviour::internals { offset += s; } return evaluators; - } // end of buildEvaluator + } // end of buildEvaluators + + static std::optional buildMassDensityEvaluator( + mgis::real& rho, MaterialStateManager& s) { + if (!isMassDensityDefined(s)) { + rho = 0; + return {}; + } + if (isMassDensityUniform(s)) { + rho = std::get(*(s.mass_density)); + return {}; + } + auto check = [&s](const auto& values) { + if (values.size() != s.n) { + mgis::raise( + "buildMassDensityEvaluator: invalid size for the arrray of the " + "mass density (" + + std::to_string(values.size()) + " values given for '" + + std::to_string(s.n) + "'integration points)"); + } + }; + if (std::holds_alternative>(*(s.mass_density))) { + const auto& values = std::get>(*(s.mass_density)); + check(values); + return std::make_tuple(0, 1, values.data()); + } + const auto& values = std::get>(*(s.mass_density)); + check(values); + return std::make_tuple(0, 1, values.data()); + } // end of buildMassDensityEvaluator static inline void applyEvaluators(std::vector& values, const std::vector& evs, @@ -147,6 +176,16 @@ namespace mgis::behaviour::internals { * end of the time step */ std::vector esvs1; + /*! + * \brief evaluator associated with the mass density at the beginning of + * the time step + */ + std::optional rho0; + /*! + * \brief evaluator associated with the mass density at the end of + * the time step + */ + std::optional rho1; }; // end of struct BehaviourEvaluators static inline BehaviourEvaluators buildBehaviourEvaluators( @@ -154,14 +193,16 @@ namespace mgis::behaviour::internals { // auto evaluators = BehaviourEvaluators{}; // treating uniform values - evaluators.mps0 = internals::buildEvaluator( + evaluators.mps0 = internals::buildEvaluators( ws.mps0, m.s0.material_properties, m, m.b.mps); - evaluators.mps1 = internals::buildEvaluator( + evaluators.mps1 = internals::buildEvaluators( ws.mps1, m.s1.material_properties, m, m.b.mps); - evaluators.esvs0 = internals::buildEvaluator( + evaluators.esvs0 = internals::buildEvaluators( ws.esvs0, m.s0.external_state_variables, m, m.b.esvs); - evaluators.esvs1 = internals::buildEvaluator( + evaluators.esvs1 = internals::buildEvaluators( ws.esvs1, m.s1.external_state_variables, m, m.b.esvs); + evaluators.rho0 = internals::buildMassDensityEvaluator(ws.rho0, m.s0); + evaluators.rho1 = internals::buildMassDensityEvaluator(ws.rho1, m.s1); return evaluators; } // end of buildBehaviourEvaluator @@ -173,6 +214,12 @@ namespace mgis::behaviour::internals { applyEvaluators(ws.mps1, evaluators.mps1, i); applyEvaluators(ws.esvs0, evaluators.esvs0, i); applyEvaluators(ws.esvs1, evaluators.esvs1, i); + if (evaluators.rho0.has_value()) { + ws.rho0 = *(std::get<2>(*(evaluators.rho0)) + i); + } + if (evaluators.rho1.has_value()) { + ws.rho1 = *(std::get<2>(*(evaluators.rho1)) + i); + } } static inline mgis::behaviour::BehaviourDataView initializeBehaviourDataView( @@ -183,6 +230,8 @@ namespace mgis::behaviour::internals { v.s1.material_properties = ws.mps1.data(); v.s0.external_state_variables = ws.esvs0.data(); v.s1.external_state_variables = ws.esvs1.data(); + v.s0.mass_density = &(ws.rho0); + v.s1.mass_density = &(ws.rho1); v.s0.stored_energy = nullptr; v.s1.stored_energy = nullptr; v.s0.dissipated_energy = nullptr; diff --git a/src/MaterialStateManager.cxx b/src/MaterialStateManager.cxx index 0f07a1a59..bff4fdb2a 100644 --- a/src/MaterialStateManager.cxx +++ b/src/MaterialStateManager.cxx @@ -212,6 +212,35 @@ namespace mgis::behaviour { return std::holds_alternative(p->second); } // end of isMaterialPropertyUniform + void setMassDensity(MaterialStateManager& m, const real v) { + m.mass_density = v; + } // end of setMassDensity + + MGIS_EXPORT void setMassDensity( + MaterialStateManager& m, + const mgis::span& v, + const MaterialStateManager::StorageMode s) { + mgis::raise_if(static_cast(v.size()) != m.n, + "setMassDensity: invalid number of values " + "(does not match the number of integration points)"); + if (s == MaterialStateManager::LOCAL_STORAGE) { + m.mass_density = std::vector{v.begin(), v.end()}; + } else { + m.mass_density = v; + } + } // end of setMassDensity + + bool isMassDensityDefined(const MaterialStateManager& m) { + return m.mass_density.has_value(); + } // end of isMassDensityDefined + + bool isMassDensityUniform(const MaterialStateManager& m) { + if (!isMassDensityDefined(m)) { + mgis::raise("isMassDensityUniform: the mass density is undefined"); + } + return std::holds_alternative(*(m.mass_density)); + } // end of isMassDensityUniform + void setExternalStateVariable(MaterialStateManager& m, const mgis::string_view& n, const real v) { @@ -365,6 +394,14 @@ namespace mgis::behaviour { for (const auto& ev : i.external_state_variables) { update_field_holder(o.external_state_variables[ev.first], ev.second); } + if (i.mass_density.has_value()) { + if (!o.mass_density.has_value()) { + o.mass_density = mgis::real{}; + } + update_field_holder(*(o.mass_density), *(i.mass_density)); + } else { + o.mass_density.reset(); + } } // end of updateValues namespace internals { diff --git a/src/State.cxx b/src/State.cxx index 96d70e3dc..4be62f6a6 100644 --- a/src/State.cxx +++ b/src/State.cxx @@ -356,6 +356,7 @@ namespace mgis::behaviour { return &v[0]; }; // end of get_ptr StateView v; + v.mass_density = &(s.mass_density); v.gradients = get_ptr(s.gradients); v.thermodynamic_forces = get_ptr(s.thermodynamic_forces); v.material_properties = get_ptr(s.material_properties); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 0486f199f..15b36eff8 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -21,7 +21,8 @@ mfront_behaviours_check_library(BehaviourTest StandardElastoViscoPlasticityPlasticityTest11 TensorialExternalStateVariableTest InitializeFunctionTest - PostProcessingTest) + PostProcessingTest + ComputeSpeedOfSound) mfront_behaviours_check_library(ModelTest ode_rk54) @@ -138,6 +139,16 @@ target_link_libraries(PostProcessingTest target_compile_definitions(PostProcessingTest PRIVATE -DTFEL_VERSION="${TFEL_VERSION}") +add_executable(ComputeSpeedOfSoundTest + EXCLUDE_FROM_ALL ComputeSpeedOfSoundTest.cxx) +target_link_libraries(ComputeSpeedOfSoundTest + PRIVATE MFrontGenericInterface) + +add_executable(ComputeSpeedOfSoundTest2 + EXCLUDE_FROM_ALL ComputeSpeedOfSoundTest2.cxx) +target_link_libraries(ComputeSpeedOfSoundTest2 + PRIVATE MFrontGenericInterface) + add_test(NAME MFrontGenericBehaviourInterfaceTest COMMAND MFrontGenericBehaviourInterfaceTest "$" "Gurson") @@ -334,3 +345,27 @@ else((CMAKE_HOST_WIN32) AND (NOT MSYS)) set_property(TEST PostProcessingTest PROPERTY DEPENDS BehaviourTest) endif((CMAKE_HOST_WIN32) AND (NOT MSYS)) + +add_test(NAME ComputeSpeedOfSoundTest + COMMAND ComputeSpeedOfSoundTest "$") +add_dependencies(check ComputeSpeedOfSoundTest) +if((CMAKE_HOST_WIN32) AND (NOT MSYS)) + set_property(TEST ComputeSpeedOfSoundTest + PROPERTY DEPENDS BehaviourTest + PROPERTY ENVIRONMENT "PATH=$\;${MGIS_PATH_STRING}") +else((CMAKE_HOST_WIN32) AND (NOT MSYS)) + set_property(TEST ComputeSpeedOfSoundTest + PROPERTY DEPENDS BehaviourTest) +endif((CMAKE_HOST_WIN32) AND (NOT MSYS)) + +add_test(NAME ComputeSpeedOfSoundTest2 + COMMAND ComputeSpeedOfSoundTest2 "$") +add_dependencies(check ComputeSpeedOfSoundTest2) +if((CMAKE_HOST_WIN32) AND (NOT MSYS)) + set_property(TEST ComputeSpeedOfSoundTest2 + PROPERTY DEPENDS BehaviourTest + PROPERTY ENVIRONMENT "PATH=$\;${MGIS_PATH_STRING}") +else((CMAKE_HOST_WIN32) AND (NOT MSYS)) + set_property(TEST ComputeSpeedOfSoundTest2 + PROPERTY DEPENDS BehaviourTest) +endif((CMAKE_HOST_WIN32) AND (NOT MSYS)) diff --git a/tests/ComputeSpeedOfSound.mfront b/tests/ComputeSpeedOfSound.mfront new file mode 100644 index 000000000..ffd02baf6 --- /dev/null +++ b/tests/ComputeSpeedOfSound.mfront @@ -0,0 +1,21 @@ +@DSL Default; +@Behaviour ComputeSpeedOfSound; + +@ProvidesSymmetricTangentOperator; +@Integrator{ + constexpr auto id = Stensor::Id(); + constexpr auto yg = stress{150e9}; + constexpr auto nu = real{0.3}; + constexpr auto l = computeLambda(yg, nu); + constexpr auto mu = computeMu(yg, nu); + const auto e = eto + deto; + sig = l * trace(e) * id + 2 * mu * e; + if (computeTangentOperator_) { + Dt = l * Stensor4::IxI() + 2 * mu * Stensor4::Id(); + } +} + +@SpeedOfSound { + constexpr auto yg = stress{150e9}; + v_sound = sqrt(yg / rho_m0); +} \ No newline at end of file diff --git a/tests/ComputeSpeedOfSoundTest.cxx b/tests/ComputeSpeedOfSoundTest.cxx new file mode 100644 index 000000000..106d8ef66 --- /dev/null +++ b/tests/ComputeSpeedOfSoundTest.cxx @@ -0,0 +1,63 @@ +/*! + * \file tests/ComputeSpeedOfSoundTest.cxx + * \brief + * \author Thomas Helfer + * \date 06/09/2023 + * \copyright (C) Copyright Thomas Helfer 2018. + * Use, modification and distribution are subject + * to one of the following licences: + * - GNU Lesser General Public License (LGPL), Version 3.0. (See accompanying + * file LGPL-3.0.txt) + * - CECILL-C, Version 1.0 (See accompanying files + * CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt). + */ + +#include +#include +#include +#include +#include "MGIS/Behaviour/State.hxx" +#include "MGIS/Behaviour/Behaviour.hxx" +#include "MGIS/Behaviour/BehaviourData.hxx" +#include "MGIS/Behaviour/Integrate.hxx" + +int main(const int argc, const char* const* argv) { + using namespace mgis; + using namespace mgis::behaviour; + if (argc != 2) { + std::cerr << "IntegrateTest: invalid number of arguments\n"; + std::exit(-1); + } + try { + const auto b = load(argv[1], "ComputeSpeedOfSound", Hypothesis::TRIDIMENSIONAL); + auto d = BehaviourData{b}; + const auto de = 5.e-5; + d.dt = 180; + // initialize the states + setExternalStateVariable(d.s1, "Temperature", 293.15); + d.s1.mass_density = 7850; + // copy d.s1 in d.s0 + update(d); + d.s1.gradients[0] = de; + // integration + d.K[0] = 100; + d.rdt = 1; + auto v = make_view(d); + integrate(v, b); + update(d); + const auto yg = 150e9; + const auto v_ref = sqrt(yg / d.s1.mass_density); + if (std::abs(d.speed_of_sound - v_ref) > 1.e-12 * v_ref) { + std::cerr << "IntegrateTest: invalid value for the speed of sound" + << "(expected '" << v_ref << "', computed '" << d.speed_of_sound + << "')\n"; + return EXIT_FAILURE; + } + } catch (std::exception& e) { + std::cerr << e.what() << '\n'; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} + + diff --git a/tests/ComputeSpeedOfSoundTest2.cxx b/tests/ComputeSpeedOfSoundTest2.cxx new file mode 100644 index 000000000..08aeebdcc --- /dev/null +++ b/tests/ComputeSpeedOfSoundTest2.cxx @@ -0,0 +1,62 @@ +/*! + * \file ComputeSpeedOfSoundTest2.cxx + * \brief + * \author Thomas Helfer + * \date 24/08/2018 + * \copyright (C) Copyright Thomas Helfer 2018. + * Use, modification and distribution are subject + * to one of the following licences: + * - GNU Lesser General Public License (LGPL), Version 3.0. (See accompanying + * file LGPL-3.0.txt) + * - CECILL-C, Version 1.0 (See accompanying files + * CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt). + */ + +#include +#include +#include +#include +#include "MGIS/Behaviour/State.hxx" +#include "MGIS/Behaviour/Behaviour.hxx" +#include "MGIS/Behaviour/MaterialDataManager.hxx" +#include "MGIS/Behaviour/Integrate.hxx" + +int main(const int argc, const char* const* argv) { + using namespace mgis; + using namespace mgis::behaviour; + const auto yg = 150e9; + const auto rho = 7850; + if (argc != 2) { + std::cerr << "ComputeSpeedOfSoundTest: invalid number of arguments\n"; + std::exit(-1); + } + try { + const auto b = load(argv[1], "ComputeSpeedOfSound", Hypothesis::TRIDIMENSIONAL); + MaterialDataManager m{b, 100}; + const auto de = 5.e-5; + // initialize the external state variable + m.s1.external_state_variables["Temperature"] = 293.15; + m.s1.mass_density = rho; + // copy d.s1 in d.s0 + update(m); + for (size_type idx = 0; idx != m.n; ++idx) { + m.s1.gradients[idx * m.s1.gradients_stride] = de; + } + const auto dt = real(180); + auto opts = BehaviourIntegrationOptions{}; + opts.integration_type = IntegrationType::INTEGRATION_NO_TANGENT_OPERATOR; + opts.compute_speed_of_sound = true; + integrate(m, opts, dt, 0, m.n); + std::cerr.precision(14); + const auto v_ref = std::sqrt(yg/rho); + if (std::abs(m.speed_of_sound[0] - v_ref) > 1.e-12 * v_ref) { + std::cerr + << "ComputeSpeedOfSoundTest: invalid value for the speed of sound\n"; + return EXIT_FAILURE; + } + } catch (std::exception& e) { + std::cerr << e.what() << '\n'; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +}