From 9b9e285894777d9d48e668074f7f0a78be7bd2c1 Mon Sep 17 00:00:00 2001 From: Thomas Helfer Date: Thu, 25 Jan 2024 10:47:01 +0100 Subject: [PATCH 1/2] Add new publications --- docs/web/bibliography.bib | 36 ++++++++++++++++++++++++++++++++++++ docs/web/publications.md | 2 ++ 2 files changed, 38 insertions(+) diff --git a/docs/web/bibliography.bib b/docs/web/bibliography.bib index 5da275fcd..c409c083d 100644 --- a/docs/web/bibliography.bib +++ b/docs/web/bibliography.bib @@ -1,3 +1,39 @@ + +@article{simo_novel_2024, + title = {A novel approach of using existing implementations of constitutive material models in any numerical codes interfacing with {MFront}}, + volume = {2}, + url = {https://www.lyellcollection.org/doi/10.1144/geoenergy2023-018}, + doi = {10.1144/geoenergy2023-018}, + abstract = {Implementing a constitutive model is a long, tedious and error-prone process, in particular for soils where a wide variety of phenomena must be taken into account. Moreover, the implementation must satisfy the interface requirements of the targeted solver. {MFront} is a popular code generator based on C++ mostly dedicated to mechanical behaviours which provides interfaces for many academic and industrial solvers. {MFront} implementations also export metadata which considerably simplifies the behaviour integration in the solver, in particular if the {MFrontGenericInterfaceSupport} ({MGIS}) is used by this solver. While {MFront} greatly reduces the amount of work required to implement a new behaviour, existing legacy implementations are highly valuable and their re-implementation should only be considered with caution considering the trade-offs. In our experience, such a re-implementation increases the maintainability and portability, and generally the numerical performances, but requires significant development effort. In this work, we developed an alternative approach, which consists in using {MFront} as a wrapper to existing legacy implementations. The {MFront} wrapper also manages the definition of appropriate metadata and handles the transfer of the data from solver to the legacy implementation on input and output. At this stage, the approach has been used to make available all constitutive models implemented in the {UMAT} format (written in Fortran) in the {OpenGeoSys} solver which is linked to {MFront} via {MGIS}. The results of a simulation using a {UMAT}-model in {OpenGeoSys} verify the approach. The usage of {MFront} as a wrapper is also shown to have an insignificant/negligible impact on the numerical performance. The proposed approach opens the door to the establishment of a new database of constitutive material models in {MFront} where legacy implementation of existing models can be made available in all solvers interfaced with {MFront}. +Thematic collection: This article is part of the Sustainable geological disposal and containment of radioactive waste collection available at: https://www.lyellcollection.org/topic/collections/radioactive}, + pages = {geoenergy2023--018}, + number = {1}, + journaltitle = {Geoenergy}, + author = {Simo, Eric and Helfer, Thomas and Mašín, David and Nagel, Thomas and Mánica, Miguel}, + urldate = {2024-01-25}, + date = {2024-12-31}, + note = {Publisher: The Geological Society of London}, +} + + +@article{singh_irradiation_2024, + title = {Irradiation dose and temperature effects on {BCC} material with dislocation based crystal plasticity}, + volume = {198}, + issn = {0306-4549}, + url = {https://www.sciencedirect.com/science/article/pii/S0306454923006060}, + doi = {10.1016/j.anucene.2023.110287}, + abstract = {Numerical modelling of temperature and strain rate dependent plasticity of body centered cubic ({BCC}) materials is carried out based on dislocation based crystal plasticity model. In view of the strong influence of radiation effects on the response of {BCC} materials, the constitutive equations are further extended to incorporate the effect of irradiation-induced defects. Size of irradiation defects produced (mainly dislocation loops) is mostly controlled by the irradiation temperature, whereas the defect number density is governed by the irradiation dose. The present constitutive model is introduced to predict the irradiation conditions and temperature dependent plasticity for {BCC} materials. The model is used to simulate the different irradiation conditions in terms of dislocation loop size and number density. To have a quantitative effect of irradiation conditions on material toughness, the Defect Induced Apparent Temperature Shift (Δ{DIAT}) is estimated and compared with available experimental evidence, for validation.}, + pages = {110287}, + journaltitle = {Annals of Nuclear Energy}, + shortjournal = {Annals of Nuclear Energy}, + author = {Singh, Kulbir and Robertson, Christian and Bhaduri, Arun Kumar}, + urldate = {2024-01-25}, + date = {2024-04-01}, + keywords = {Crystal plasticity, Dislocation mobility, Irradiation hardening, Screw dislocations}, + file = {ScienceDirect Snapshot:/home/th202608/.zotero/zotero/vzpzva46.default/zotero/storage/KBSDHI6Z/S0306454923006060.html:text/html}, +} + + @article{senac_yield_2023, title = {Yield surface for void growth and coalescence of porous anisotropic materials under axisymmetric loading}, volume = {179}, diff --git a/docs/web/publications.md b/docs/web/publications.md index bfdb062dd..b832aa817 100644 --- a/docs/web/publications.md +++ b/docs/web/publications.md @@ -13,6 +13,8 @@ eqnPrefixTemplate: "($$i$$)" --- nocite: | + @simo_novel_2024 + @singh_irradiation_2024 @fokam_implementation_2023 @bacquaert_standard_2023 @rapanakis_three-dimensional_2023 From af4c414ba60a23fdf76a49bfb740fa0aac4fbc6b Mon Sep 17 00:00:00 2001 From: Thomas Helfer Date: Thu, 25 Jan 2024 11:37:04 +0100 Subject: [PATCH 2/2] Introduce the ST2toST2Concept. Remove usage of implementsST2toST2Concept. Code clean-up --- include/CMakeLists.txt | 1 - .../FiniteStrainBehaviourTangentOperator.hxx | 7 +- include/TFEL/Math/ST2toST2/ChangeBasis.hxx | 97 --- ...onvertLogarithmicStrainTangentOperator.hxx | 37 +- ...atialModuliToKirchhoffJaumanRateModuli.hxx | 199 +++--- .../ST2toST2/ConvertT2toST2ToST2toST2Expr.hxx | 6 +- .../TFEL/Math/ST2toST2/ST2toST2Concept.hxx | 174 ++---- .../TFEL/Math/ST2toST2/ST2toST2Concept.ixx | 103 ++-- .../TFEL/Math/ST2toST2/ST2toST2ConceptIO.hxx | 6 +- .../ST2toST2/ST2toST2ConceptOperations.hxx | 10 +- .../ST2toST2/ST2toST2ConceptPushForward.ixx | 566 +++++++++--------- .../ST2toST2/ST2toST2ST2toST2ProductExpr.hxx | 18 +- .../ST2toST2/ST2toST2StensorProductExpr.hxx | 61 +- .../Math/ST2toST2/ST2toST2TransposeExpr.hxx | 26 +- .../ST2toST2/StensorST2toST2ProductExpr.hxx | 84 +-- .../Math/ST2toST2/StensorSquareDerivative.hxx | 70 +-- include/TFEL/Math/ST2toST2/st2tost2.ixx | 273 +++++---- .../ST2toT2/ST2toT2ST2toST2ProductExpr.hxx | 36 +- .../StensorProductLeftDerivativeExpr.hxx | 50 +- .../StensorProductRightDerivativeExpr.hxx | 44 +- .../ST2toT2/T2toST2ST2toT2ProductExpr.hxx | 6 +- include/TFEL/Math/ST2toT2/st2tot2.ixx | 28 +- ...ecompositionInPositiveAndNegativeParts.hxx | 33 +- ...ecompositionInPositiveAndNegativeParts.ixx | 33 +- .../Stensor/Internals/StensorChangeBasis.hxx | 11 +- .../StensorComputeEigenVectorsDerivatives.hxx | 63 +- ...nsorComputeIsotropicFunctionDerivative.hxx | 123 ++-- include/TFEL/Math/Stensor/StensorConcept.hxx | 4 +- include/TFEL/Math/Stensor/stensor.ixx | 23 +- .../T2toST2/ST2toST2T2toST2ProductExpr.hxx | 62 +- include/TFEL/Math/Tensor/tensor.ixx | 26 +- include/TFEL/Math/st2tost2.hxx | 146 ++--- include/TFEL/Math/st2tot2.hxx | 21 +- include/TFEL/Math/stensor.hxx | 24 +- include/TFEL/Math/tensor.hxx | 20 +- 35 files changed, 1099 insertions(+), 1392 deletions(-) delete mode 100644 include/TFEL/Math/ST2toST2/ChangeBasis.hxx diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt index 3e02b575b..abdfcb937 100644 --- a/include/CMakeLists.txt +++ b/include/CMakeLists.txt @@ -397,7 +397,6 @@ install_header(TFEL/Math/ST2toST2 ST2toST2ST2toST2ProductExpr.hxx) install_header(TFEL/Math/ST2toST2 ConvertT2toST2ToST2toST2Expr.hxx) install_header(TFEL/Math/ST2toST2 StensorSquareDerivative.hxx) install_header(TFEL/Math/ST2toST2 BuildFromRotationMatrix.hxx) -install_header(TFEL/Math/ST2toST2 ChangeBasis.hxx) install_header(TFEL/Math/ST2toST2 StensorSymmetricProductDerivative.hxx) install_header(TFEL/Math/ST2toST2 ConvertToTangentModuli.hxx) install_header(TFEL/Math/ST2toST2 ConvertSpatialModuliToKirchhoffJaumanRateModuli.hxx) diff --git a/include/TFEL/Material/FiniteStrainBehaviourTangentOperator.hxx b/include/TFEL/Material/FiniteStrainBehaviourTangentOperator.hxx index f8e25ac51..fc7fafe65 100644 --- a/include/TFEL/Material/FiniteStrainBehaviourTangentOperator.hxx +++ b/include/TFEL/Material/FiniteStrainBehaviourTangentOperator.hxx @@ -149,12 +149,11 @@ namespace tfel::material { return *this; } /*! - * assignement operator + * \brief assignement operator */ - template + template std::enable_if_t< - tfel::math::implementsST2toST2Concept() && - tfel::math::getSpaceDimension() == N && + tfel::math::getSpaceDimension() == N && std::is_same_v, StressType>, FiniteStrainBehaviourTangentOperator&> operator=(const T& e) { diff --git a/include/TFEL/Math/ST2toST2/ChangeBasis.hxx b/include/TFEL/Math/ST2toST2/ChangeBasis.hxx deleted file mode 100644 index d4f000072..000000000 --- a/include/TFEL/Math/ST2toST2/ChangeBasis.hxx +++ /dev/null @@ -1,97 +0,0 @@ -/*! - * \file include/TFEL/Math/ST2toST2/ChangeBasis.hxx - * \brief - * \author Thomas Helfer - * \date 13 oct. 2015 - * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights - * reserved. - * This project is publicly released under either the GNU GPL Licence - * or the CECILL-A licence. A copy of thoses licences are delivered - * with the sources of TFEL. CEA or EDF may also distribute this - * project under specific licensing conditions. - */ - -#ifndef LIB_TFEL_MATH_ST2TOST2_CHANGEBASIS_HXX -#define LIB_TFEL_MATH_ST2TOST2_CHANGEBASIS_HXX - -#include - -namespace tfel::math::st2tost2_internals { - - /*! - * \brief an auxiliary structure used to dispatch over the - * dimension - */ - template - struct ChangeBasis; - - /*! - * \brief partial specialisation in 1D - */ - template <> - struct ChangeBasis<1u> { - /*! - */ - template - static TFEL_MATH_INLINE2 typename std::enable_if< - (tfel::math::implementsST2toST2Concept() && - tfel::math::getSpaceDimension() == 1u), - tfel::math::st2tost2<1u, numeric_type>>::type - exe(const ST2toST2Type& s, - const tfel::math::rotation_matrix< - tfel::math::numeric_type>&) { - return s; - } - }; // end of ChangeBasis<1u> - /*! - * \brief partial specialisation in 2D - */ - template <> - struct ChangeBasis<2u> { - /*! - * \param[in] r : rotation matrix - */ - template - static TFEL_MATH_INLINE2 typename std::enable_if< - (tfel::math::implementsST2toST2Concept() && - tfel::math::getSpaceDimension() == 2u), - tfel::math::st2tost2<2u, tfel::math::numeric_type>>::type - exe(const ST2toST2Type& s, - const tfel::math::rotation_matrix< - tfel::math::numeric_type>& r) { - using st2tost2 = - tfel::math::st2tost2<2u, tfel::math::numeric_type>; - const auto sr = st2tost2::fromRotationMatrix(r); - const auto sir = st2tost2::fromRotationMatrix(transpose(r)); - return sr * s * sir; - } // end of ChangeBasis<2u>::exe - }; // end of ChangeBasis<2u> - /*! - * \brief partial specialisation in 3D - */ - template <> - struct ChangeBasis<3u> { - /*! - * \return the st2tost2 that has the same effect as applying - * the rotation - * \param[in] r : rotation matrix - */ - template - static TFEL_MATH_INLINE2 typename std::enable_if< - (tfel::math::implementsST2toST2Concept() && - tfel::math::getSpaceDimension() == 3u), - tfel::math::st2tost2<3u, tfel::math::numeric_type>>::type - exe(const ST2toST2Type& s, - const tfel::math::rotation_matrix< - tfel::math::numeric_type>& r) { - using st2tost2 = - tfel::math::st2tost2<3u, tfel::math::numeric_type>; - const auto sr = st2tost2::fromRotationMatrix(r); - const auto sir = st2tost2::fromRotationMatrix(transpose(r)); - return sr * s * sir; - } // end of ChangeBasis<3u>::exe - }; // end of ChangeBasis<2u> - -} // end of namespace tfel::math::st2tost2_internals - -#endif /* LIB_TFEL_MATH_ST2TOST2_CHANGEBASIS_HXX */ diff --git a/include/TFEL/Math/ST2toST2/ConvertLogarithmicStrainTangentOperator.hxx b/include/TFEL/Math/ST2toST2/ConvertLogarithmicStrainTangentOperator.hxx index 7b37b23f9..45309d4f8 100644 --- a/include/TFEL/Math/ST2toST2/ConvertLogarithmicStrainTangentOperator.hxx +++ b/include/TFEL/Math/ST2toST2/ConvertLogarithmicStrainTangentOperator.hxx @@ -54,27 +54,24 @@ namespace tfel::math { * \param[in] m: eigen vectors of the right Cauchy Green tensor * \param[in] vp: eigen values of the right Cauchy Green tensor */ - template - static typename std::enable_if< - ((implementsST2toST2Concept()) && - (tfel::math::getSpaceDimension() == 1u) && - (std::is_same, stress>::value) && - (implementsST2toST2Concept()) && - (tfel::math::getSpaceDimension() == 1u) && - (isAssignableTo, stress>()), - (implementsST2toST2Concept()) && - (tfel::math::getSpaceDimension() == 1u) && - (isAssignableTo, real>())), - void>::type - exe(ST2toST2Type& Cse, - const ST2toST2Type2& C, - const ST2toST2Type2& P, - const StressStensorType& T, - const tmatrix<3u, 3u, real>&, - const tvector<3u, real>&) { + TFEL_HOST_DEVICE constexpr void exe(ST2toST2Type& Cse, + const ST2toST2Type2& C, + const ST2toST2Type2& P, + const StressStensorType& T, + const tmatrix<3u, 3u, real>&, + const tvector<3u, + real>&) noexcept // + requires( + (tfel::math::getSpaceDimension() == 1u) && + (std::is_same, stress>::value) && + (tfel::math::getSpaceDimension() == 1u) && + (isAssignableTo, stress>()), + (tfel::math::getSpaceDimension() == 1u) && + (isAssignableTo, real>())) { const auto iC0 = P(0, 0); const auto iC1 = P(1, 1); const auto iC2 = P(2, 2); diff --git a/include/TFEL/Math/ST2toST2/ConvertSpatialModuliToKirchhoffJaumanRateModuli.hxx b/include/TFEL/Math/ST2toST2/ConvertSpatialModuliToKirchhoffJaumanRateModuli.hxx index f48ceb602..3680a9341 100644 --- a/include/TFEL/Math/ST2toST2/ConvertSpatialModuliToKirchhoffJaumanRateModuli.hxx +++ b/include/TFEL/Math/ST2toST2/ConvertSpatialModuliToKirchhoffJaumanRateModuli.hxx @@ -32,124 +32,89 @@ namespace tfel::math { * \return the moduli associated with Jauman rate of the Kirchhoff * stress */ - template - typename std::enable_if< - ((getSpaceDimension() == 1u) && - (getSpaceDimension() == 1u) && - implementsST2toST2Concept() && - std::is_same_v, numeric_type>), - st2tost2<1u, numeric_type>>::type - convertSpatialModuliToKirchhoffJaumanRateModuli(const ST2toST2Type& C_s, - const StensorType& tau) { + template + TFEL_HOST_DEVICE constexpr auto + convertSpatialModuliToKirchhoffJaumanRateModuli( + const ST2toST2Type& C_s, + const StensorType& + tau) noexcept requires((getSpaceDimension() == + getSpaceDimension()) && + (std::is_same_v, + numeric_type>)) { + constexpr auto N = getSpaceDimension(); using NumType = numeric_type; - st2tost2<1u, NumType> C_tJ; - C_tJ(0, 0) = C_s(0, 0) + 2 * tau[0]; - C_tJ(0, 1) = C_s(0, 1); - C_tJ(0, 2) = C_s(0, 2); - C_tJ(1, 0) = C_s(1, 0); - C_tJ(1, 1) = C_s(1, 1) + 2 * tau[1]; - C_tJ(1, 2) = C_s(1, 2); - C_tJ(2, 0) = C_s(2, 0); - C_tJ(2, 1) = C_s(2, 1); - C_tJ(2, 2) = C_s(2, 2) + 2 * tau[2]; - return C_tJ; - } - - /*! - * \brief convert the spatial moduli in the moduli associated with - * Jauman rate of the Kirchhoff stress - * \param[in] C_s: spatial moduli - * \param[in] tau: Kirchhoff stress - * \return the moduli associated with Jauman rate of the Kirchhoff - * stress - */ - template - typename std::enable_if< - ((getSpaceDimension() == 2u) && - (getSpaceDimension() == 2u) && - implementsST2toST2Concept() && - std::is_same_v, numeric_type>), - st2tost2<2u, numeric_type>>::type - convertSpatialModuliToKirchhoffJaumanRateModuli(const ST2toST2Type& C_s, - const StensorType& tau) { - using NumType = numeric_type; - st2tost2<2u, NumType> C_tJ; - C_tJ(0, 0) = C_s(0, 0) + 2 * tau[0]; - C_tJ(1, 1) = C_s(1, 1) + 2 * tau[1]; - C_tJ(2, 2) = C_s(2, 2) + 2 * tau[2]; - C_tJ(3, 0) = C_s(3, 0) + tau[3]; - C_tJ(3, 1) = C_s(3, 1) + tau[3]; - C_tJ(0, 3) = C_s(0, 3) + tau[3]; - C_tJ(1, 3) = C_s(1, 3) + tau[3]; - C_tJ(3, 3) = C_s(3, 3) + tau[1] + tau[0]; - C_tJ(0, 1) = C_s(0, 1); - C_tJ(0, 2) = C_s(0, 2); - C_tJ(3, 2) = C_s(3, 2); - C_tJ(1, 0) = C_s(1, 0); - C_tJ(1, 2) = C_s(1, 2); - C_tJ(2, 0) = C_s(2, 0); - C_tJ(2, 1) = C_s(2, 1); - C_tJ(2, 3) = C_s(2, 3); - return C_tJ; - } - - /*! - * \brief convert the spatial moduli in the moduli associated with - * Jauman rate of the Kirchhoff stress - * \param[in] C_s: spatial moduli - * \param[in] tau: Kirchhoff stress - * \return the moduli associated with Jauman rate of the Kirchhoff - * stress - */ - template - typename std::enable_if< - ((getSpaceDimension() == 3u) && - (getSpaceDimension() == 3u) && - implementsST2toST2Concept() && - std::is_same_v, numeric_type>), - st2tost2<3u, numeric_type>>::type - convertSpatialModuliToKirchhoffJaumanRateModuli(const ST2toST2Type& C_s, - const StensorType& tau) { - using NumType = numeric_type; - constexpr auto icste = Cste::isqrt2; - st2tost2<3u, NumType> C_tJ; - C_tJ(0, 0) = C_s(0, 0) + 2 * tau[0]; - C_tJ(1, 1) = C_s(1, 1) + 2 * tau[1]; - C_tJ(2, 2) = C_s(2, 2) + 2 * tau[2]; - C_tJ(0, 3) = C_s(0, 3) + tau[3]; - C_tJ(0, 4) = C_s(0, 4) + tau[4]; - C_tJ(1, 3) = C_s(1, 3) + tau[3]; - C_tJ(1, 5) = C_s(1, 5) + tau[5]; - C_tJ(2, 4) = C_s(2, 4) + tau[4]; - C_tJ(2, 5) = C_s(2, 5) + tau[5]; - C_tJ(3, 3) = C_s(3, 3) + tau[1] + tau[0]; - C_tJ(3, 0) = C_s(3, 0) + tau[3]; - C_tJ(3, 1) = C_s(3, 1) + tau[3]; - C_tJ(3, 4) = C_s(3, 4) + tau[5] * icste; - C_tJ(3, 5) = C_s(3, 5) + tau[4] * icste; - C_tJ(4, 0) = C_s(4, 0) + tau[4]; - C_tJ(4, 2) = C_s(4, 2) + tau[4]; - C_tJ(4, 4) = C_s(4, 4) + tau[2] + tau[0]; - C_tJ(4, 3) = C_s(4, 3) + tau[5] * icste; - C_tJ(4, 5) = C_s(4, 5) + tau[3] * icste; - C_tJ(5, 1) = C_s(5, 1) + tau[5]; - C_tJ(5, 2) = C_s(5, 2) + tau[5]; - C_tJ(5, 3) = C_s(5, 3) + tau[4] * icste; - C_tJ(5, 4) = C_s(5, 4) + tau[3] * icste; - C_tJ(5, 5) = C_s(5, 5) + tau[2] + tau[1]; - C_tJ(0, 1) = C_s(0, 1); - C_tJ(0, 2) = C_s(0, 2); - C_tJ(5, 0) = C_s(5, 0); - C_tJ(0, 5) = C_s(0, 5); - C_tJ(1, 0) = C_s(1, 0); - C_tJ(1, 2) = C_s(1, 2); - C_tJ(1, 4) = C_s(1, 4); - C_tJ(2, 0) = C_s(2, 0); - C_tJ(2, 1) = C_s(2, 1); - C_tJ(2, 3) = C_s(2, 3); - C_tJ(3, 2) = C_s(3, 2); - C_tJ(4, 1) = C_s(4, 1); - return C_tJ; + if constexpr (N == 1) { + st2tost2<1u, NumType> C_tJ; + C_tJ(0, 0) = C_s(0, 0) + 2 * tau[0]; + C_tJ(0, 1) = C_s(0, 1); + C_tJ(0, 2) = C_s(0, 2); + C_tJ(1, 0) = C_s(1, 0); + C_tJ(1, 1) = C_s(1, 1) + 2 * tau[1]; + C_tJ(1, 2) = C_s(1, 2); + C_tJ(2, 0) = C_s(2, 0); + C_tJ(2, 1) = C_s(2, 1); + C_tJ(2, 2) = C_s(2, 2) + 2 * tau[2]; + return C_tJ; + } else if constexpr (N == 2) { + st2tost2<2u, NumType> C_tJ; + C_tJ(0, 0) = C_s(0, 0) + 2 * tau[0]; + C_tJ(1, 1) = C_s(1, 1) + 2 * tau[1]; + C_tJ(2, 2) = C_s(2, 2) + 2 * tau[2]; + C_tJ(3, 0) = C_s(3, 0) + tau[3]; + C_tJ(3, 1) = C_s(3, 1) + tau[3]; + C_tJ(0, 3) = C_s(0, 3) + tau[3]; + C_tJ(1, 3) = C_s(1, 3) + tau[3]; + C_tJ(3, 3) = C_s(3, 3) + tau[1] + tau[0]; + C_tJ(0, 1) = C_s(0, 1); + C_tJ(0, 2) = C_s(0, 2); + C_tJ(3, 2) = C_s(3, 2); + C_tJ(1, 0) = C_s(1, 0); + C_tJ(1, 2) = C_s(1, 2); + C_tJ(2, 0) = C_s(2, 0); + C_tJ(2, 1) = C_s(2, 1); + C_tJ(2, 3) = C_s(2, 3); + return C_tJ; + } else { + constexpr auto icste = Cste::isqrt2; + st2tost2<3u, NumType> C_tJ; + C_tJ(0, 0) = C_s(0, 0) + 2 * tau[0]; + C_tJ(1, 1) = C_s(1, 1) + 2 * tau[1]; + C_tJ(2, 2) = C_s(2, 2) + 2 * tau[2]; + C_tJ(0, 3) = C_s(0, 3) + tau[3]; + C_tJ(0, 4) = C_s(0, 4) + tau[4]; + C_tJ(1, 3) = C_s(1, 3) + tau[3]; + C_tJ(1, 5) = C_s(1, 5) + tau[5]; + C_tJ(2, 4) = C_s(2, 4) + tau[4]; + C_tJ(2, 5) = C_s(2, 5) + tau[5]; + C_tJ(3, 3) = C_s(3, 3) + tau[1] + tau[0]; + C_tJ(3, 0) = C_s(3, 0) + tau[3]; + C_tJ(3, 1) = C_s(3, 1) + tau[3]; + C_tJ(3, 4) = C_s(3, 4) + tau[5] * icste; + C_tJ(3, 5) = C_s(3, 5) + tau[4] * icste; + C_tJ(4, 0) = C_s(4, 0) + tau[4]; + C_tJ(4, 2) = C_s(4, 2) + tau[4]; + C_tJ(4, 4) = C_s(4, 4) + tau[2] + tau[0]; + C_tJ(4, 3) = C_s(4, 3) + tau[5] * icste; + C_tJ(4, 5) = C_s(4, 5) + tau[3] * icste; + C_tJ(5, 1) = C_s(5, 1) + tau[5]; + C_tJ(5, 2) = C_s(5, 2) + tau[5]; + C_tJ(5, 3) = C_s(5, 3) + tau[4] * icste; + C_tJ(5, 4) = C_s(5, 4) + tau[3] * icste; + C_tJ(5, 5) = C_s(5, 5) + tau[2] + tau[1]; + C_tJ(0, 1) = C_s(0, 1); + C_tJ(0, 2) = C_s(0, 2); + C_tJ(5, 0) = C_s(5, 0); + C_tJ(0, 5) = C_s(0, 5); + C_tJ(1, 0) = C_s(1, 0); + C_tJ(1, 2) = C_s(1, 2); + C_tJ(1, 4) = C_s(1, 4); + C_tJ(2, 0) = C_s(2, 0); + C_tJ(2, 1) = C_s(2, 1); + C_tJ(2, 3) = C_s(2, 3); + C_tJ(3, 2) = C_s(3, 2); + C_tJ(4, 1) = C_s(4, 1); + return C_tJ; + } } } // end of namespace tfel::math diff --git a/include/TFEL/Math/ST2toST2/ConvertT2toST2ToST2toST2Expr.hxx b/include/TFEL/Math/ST2toST2/ConvertT2toST2ToST2toST2Expr.hxx index 348c885a4..0a4110636 100644 --- a/include/TFEL/Math/ST2toST2/ConvertT2toST2ToST2toST2Expr.hxx +++ b/include/TFEL/Math/ST2toST2/ConvertT2toST2ToST2toST2Expr.hxx @@ -31,7 +31,7 @@ namespace tfel::math { */ template struct Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder<9u, numeric_type> { static_assert(getSpaceDimension() == 1u); @@ -77,7 +77,7 @@ namespace tfel::math { */ template struct Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder<16u, numeric_type> { //! a simple check @@ -140,7 +140,7 @@ namespace tfel::math { */ template struct Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder<36u, numeric_type> { static_assert(getSpaceDimension() == 3u); diff --git a/include/TFEL/Math/ST2toST2/ST2toST2Concept.hxx b/include/TFEL/Math/ST2toST2/ST2toST2Concept.hxx index f2346bbaa..7a86c2f68 100644 --- a/include/TFEL/Math/ST2toST2/ST2toST2Concept.hxx +++ b/include/TFEL/Math/ST2toST2/ST2toST2Concept.hxx @@ -27,48 +27,61 @@ namespace tfel::math { - //! a simple alias - template - struct ST2toST2TransposeExpr; - /*! * \class ST2toST2Tag * \brief Helper class to characterise st2tost2. */ struct ST2toST2Tag {}; - + /*! + * \brief an helper class that simply exposes publically a member named + * ConceptTag as an alias to StensorTag. + * + * The main reason for this alias is to properly implement the `ConceptRebind` + * metafunction. + */ template - struct ST2toST2Concept { - typedef ST2toST2Tag ConceptTag; - - protected: - ST2toST2Concept() = default; - ST2toST2Concept(ST2toST2Concept&&) = default; - ST2toST2Concept(const ST2toST2Concept&) = default; - ST2toST2Concept& operator=(const ST2toST2Concept&) = default; - ~ST2toST2Concept() = default; + struct ST2toST2ConceptBase { + using ConceptTag = ST2toST2Tag; }; - /*! - * \brief an helper function which returns if the given type implements the - * `ST2toST2Concept`. - * \tparam ST2toST2Type: type tested + * \brief definition of the StensorConcept + * a class matching the stensor concept must expose the `StensorTag` and have + * access operators. */ - template - TFEL_HOST_DEVICE constexpr bool implementsST2toST2Concept() { - return tfel::meta::implements(); - } // end of implementsST2toST2Concept + template + concept ST2toST2Concept = + (std::is_same_v::ConceptTag, ST2toST2Tag>)&& // + (requires(const T t, const unsigned short i, const unsigned short j) { + t(i, j); + }); + //! a simple alias + template + struct ST2toST2TransposeExpr; + //! \brief partial specialisation for symmetric tensors template struct ConceptRebind { - typedef ST2toST2Concept type; + //! \brief a simple alias + using type = ST2toST2ConceptBase; }; - template - std::enable_if_t< - implementsST2toST2Concept(), - typename tfel::typetraits::AbsType>::type> - abs(const ST2toST2Type&); + /*! + * \return the sum of the absolute values of all components of an + * linear application transforming a symmetric tensor in a symmetric tensor + * \param[in] s: linear application transforming a symmetric tensor in a + * symmetric tensor + */ + TFEL_HOST_DEVICE constexpr auto abs(const ST2toST2Concept auto&) noexcept; + /*! + * \return a transposed view of a fourth order tensor + * \param[in] t: fourth order tensor + */ + TFEL_HOST_DEVICE constexpr auto transpose(ST2toST2Concept auto&&) noexcept; + /*! + * \return the determinant of a `st2tost2` + * \param[in] s: fourth order tensor + */ + TFEL_HOST_DEVICE constexpr auto det(const ST2toST2Concept auto&) noexcept; /*! * \brief compute de derivative of the push-forward of a symmetric @@ -79,12 +92,12 @@ namespace tfel::math { * \param[out] r: derivative of the push-forward symmetric tensor * \param[in] F: deformation gradient */ - template - typename std::enable_if() && - tfel::typetraits::IsFundamentalNumericType< - numeric_type>::cond, - void>::type - computePushForwardDerivative(ST2toST2ResultType&, const TensorType&); + template + TFEL_HOST_DEVICE constexpr void computePushForwardDerivative( + ST2toST2ResultType&, + const TensorType&) noexcept // + requires(tfel::typetraits::IsFundamentalNumericType< + numeric_type>::cond); /*! * \brief performs the push_forward of a st2tost2: * \[ @@ -94,85 +107,28 @@ namespace tfel::math { * \param[in] C: input * \param[in] F: deformation gradient */ - template - std::enable_if_t() && - implementsST2toST2Concept() && - getSpaceDimension() == 1u && - getSpaceDimension() == 1u && - getSpaceDimension() == 1u, - void> - push_forward(ST2toST2Type&, const ST2toST2Type2&, const TensorType&); - /*! - * \brief performs the push_forward of a st2tost2: - * \[ - * Ct_{ijkl}=F_{im}F_{jn}F_{kp}F_{lq}C_{mnpq} - * \] - * \param[out] Ct: result - * \param[in] C: input - * \param[in] F: deformation gradient - */ - template - std::enable_if_t() && - implementsST2toST2Concept() && - getSpaceDimension() == 2u && - getSpaceDimension() == 2u && - getSpaceDimension() == 2u, - void> - push_forward(ST2toST2Type&, const ST2toST2Type2&, const TensorType&); - /*! - * \brief performs the push_forward of a st2tost2: - * \[ - * Ct_{ijkl}=F_{im}F_{jn}F_{kp}F_{lq}C_{mnpq} - * \] - * \param[out] Ct: result - * \param[in] C: input - * \param[in] F: deformation gradient - */ - template - std::enable_if_t() && - implementsST2toST2Concept() && - getSpaceDimension() == 3u && - getSpaceDimension() == 3u && - getSpaceDimension() == 3u, - void> - push_forward(ST2toST2Type&, const ST2toST2Type2&, const TensorType&); - /*! - * \return a transposed view - */ - template - TFEL_MATH_INLINE auto transpose(ST2toST2Type&& t) -> std::enable_if_t< - implementsST2toST2Concept(), - Expr, ST2toST2TransposeExpr>>; - /*! - * \return the determinant of a `st2tost2` - * \param[in] s: fourth order tensor - */ - template - std::enable_if_t< - implementsST2toST2Concept() && - (getSpaceDimension() == 1u) && - isScalar>(), - typename ComputeUnaryResult, Power<3>>::Result> - det(const ST2toST2Type&); + TFEL_HOST_DEVICE constexpr void push_forward(ST2toST2Type&, + const ST2toST2Type2&, + const TensorType&) noexcept // + requires(getSpaceDimension() == + getSpaceDimension() && + getSpaceDimension() == + getSpaceDimension()); + /*! - * \return the determinant of a `st2tost2` - * \param[in] s: fourth order tensor + * \brief an helper function which returns if the given type implements the + * `ST2toST2Concept`. + * \tparam ST2toST2Type: type tested + * \note function given for backward compatibility with versions prior + * to 5.0 */ template - std::enable_if_t() && - ((getSpaceDimension() == 2u) || - (getSpaceDimension() == 3u)) && - isScalar>(), - typename ComputeUnaryResult< - numeric_type, - Power()>>::Result> - det(const ST2toST2Type&); + [[deprecated]] TFEL_HOST_DEVICE constexpr bool implementsST2toST2Concept() { + return ST2toST2Concept; + } // end of implementsST2toST2Concept } // end of namespace tfel::math diff --git a/include/TFEL/Math/ST2toST2/ST2toST2Concept.ixx b/include/TFEL/Math/ST2toST2/ST2toST2Concept.ixx index 26ba0b41f..13d3d1995 100644 --- a/include/TFEL/Math/ST2toST2/ST2toST2Concept.ixx +++ b/include/TFEL/Math/ST2toST2/ST2toST2Concept.ixx @@ -21,86 +21,69 @@ namespace tfel::math { - template - std::enable_if_t< - implementsST2toST2Concept(), - typename tfel::typetraits::AbsType>::type> - abs(const ST2toST2Type& v) { + TFEL_HOST_DEVICE constexpr auto abs(const ST2toST2Concept auto& v) noexcept { + using ST2toST2Type = decltype(v); using NumType = numeric_type; using IndexType = index_type; using AbsNumType = typename tfel::typetraits::AbsType::type; constexpr auto ssize = StensorDimeToSize()>::value; - AbsNumType a(0); + auto a = AbsNumType{}; for (IndexType i = 0; i < ssize; ++i) { for (IndexType j = 0; j < ssize; ++j) { a += abs(v(i, j)); } } return a; - } + } // end of abs - template - auto transpose(ST2toST2Type&& t) - -> std::enable_if_t(), - Expr, - ST2toST2TransposeExpr>> { + TFEL_HOST_DEVICE constexpr auto transpose(ST2toST2Concept auto&& t) noexcept { + using ST2toST2Type = decltype(t); return Expr, ST2toST2TransposeExpr>( std::forward(t)); - } + } // end of transpose - template - std::enable_if_t< - implementsST2toST2Concept() && - (getSpaceDimension() == 1u) && - isScalar>(), - typename ComputeUnaryResult, Power<3>>::Result> - det(const ST2toST2Type& s) { - const auto a = s(0, 0); - const auto b = s(0, 1); - const auto c = s(0, 2); - const auto d = s(1, 0); - const auto e = s(1, 1); - const auto f = s(1, 2); - const auto g = s(2, 0); - const auto h = s(2, 1); - const auto i = s(2, 2); - return a * (e * i - f * h) + b * (f * g - d * i) + c * (d * h - e * g); - } // end of det - - template - std::enable_if_t() && - ((getSpaceDimension() == 2u) || - (getSpaceDimension() == 3u)) && - isScalar>(), - typename ComputeUnaryResult< - numeric_type, - Power()>>::Result> - det(const ST2toST2Type& s) { - using real = numeric_type; + TFEL_HOST_DEVICE constexpr auto det(const ST2toST2Concept auto& s) noexcept { + using ST2toST2Type = decltype(s); constexpr auto N = getSpaceDimension(); - constexpr auto ts = StensorDimeToSize::value; - tmatrix m; - TinyPermutation p; - tfel::fsalgo::copy::exe(s.begin(), m.begin()); - const auto r = LUDecomp::exe(m, p); - if (!r.first) { - return {}; - } - auto v = base_type{1}; - for (const index_type i = 0; i != ts; ++i) { - v *= m(i, i); + if constexpr (N == 1) { + const auto a = s(0, 0); + const auto b = s(0, 1); + const auto c = s(0, 2); + const auto d = s(1, 0); + const auto e = s(1, 1); + const auto f = s(1, 2); + const auto g = s(2, 0); + const auto h = s(2, 1); + const auto i = s(2, 2); + return a * (e * i - f * h) + b * (f * g - d * i) + c * (d * h - e * g); + } else { + constexpr auto ts = StensorDimeToSize::value; + using Result = UnaryResultType, Power>; + using real = base_type>; + tmatrix m; + tfel::fsalgo::transform::exe( + s.begin(), m.begin(), [](const auto v) { return base_type_cast(v); }); + TinyPermutation p; + const auto r = LUDecomp::exe(m, p); + if (!r.first) { + return Result{}; + } + auto v = base_type{1}; + for (const index_type i = 0; i != ts; ++i) { + v *= m(i, i); + } + return r.second == 1 ? Result{v} : -Result{v}; } - return r.second == 1 ? v : -v; } // end of det - template - typename std::enable_if() && - tfel::typetraits::IsFundamentalNumericType< - numeric_type>::cond, - void>::type - computePushForwardDerivative(ST2toST2ResultType& r, const TensorType& F) { + template + TFEL_HOST_DEVICE constexpr void computePushForwardDerivative( + ST2toST2ResultType& r, + const TensorType& F) noexcept // + requires(tfel::typetraits::IsFundamentalNumericType< + numeric_type>::cond) { constexpr auto N = getSpaceDimension(); static_assert(getSpaceDimension() == N); using value_type = numeric_type; diff --git a/include/TFEL/Math/ST2toST2/ST2toST2ConceptIO.hxx b/include/TFEL/Math/ST2toST2/ST2toST2ConceptIO.hxx index 8fbd458e5..2163565f2 100644 --- a/include/TFEL/Math/ST2toST2/ST2toST2ConceptIO.hxx +++ b/include/TFEL/Math/ST2toST2/ST2toST2ConceptIO.hxx @@ -20,9 +20,9 @@ namespace tfel::math { // Serialisation operator - template - std::enable_if_t(), std::ostream&> - operator<<(std::ostream& os, const ST2toST2Type& s) { + TFEL_HOST std::ostream& operator<<(std::ostream& os, + const ST2toST2Concept auto& s) { + using ST2toST2Type = decltype(s); constexpr auto stensor_size = StensorDimeToSize()>::value; os << "["; diff --git a/include/TFEL/Math/ST2toST2/ST2toST2ConceptOperations.hxx b/include/TFEL/Math/ST2toST2/ST2toST2ConceptOperations.hxx index 2783b40c8..68130e91e 100644 --- a/include/TFEL/Math/ST2toST2/ST2toST2ConceptOperations.hxx +++ b/include/TFEL/Math/ST2toST2/ST2toST2ConceptOperations.hxx @@ -97,15 +97,11 @@ namespace tfel::math { Expr>>; }; - template - TFEL_MATH_INLINE typename std::enable_if< - implementsST2toST2Concept() && - !isInvalid>(), - BinaryOperationHandler>::type - operator|(const T1& a, const T2& b) { + template + TFEL_HOST_DEVICE constexpr auto operator|(const T1& a, const T2& b) noexcept { using Handle = BinaryOperationHandler; return Handle(a, b); - } + } // end of operator| } // end of namespace tfel::math diff --git a/include/TFEL/Math/ST2toST2/ST2toST2ConceptPushForward.ixx b/include/TFEL/Math/ST2toST2/ST2toST2ConceptPushForward.ixx index 4c76d4fc7..a6e6a13d0 100644 --- a/include/TFEL/Math/ST2toST2/ST2toST2ConceptPushForward.ixx +++ b/include/TFEL/Math/ST2toST2/ST2toST2ConceptPushForward.ixx @@ -18,312 +18,292 @@ namespace tfel::math { - template - typename std::enable_if() && - implementsST2toST2Concept() && - getSpaceDimension() == 1u && - getSpaceDimension() == 1u && - getSpaceDimension() == 1u, - void>::type - push_forward(ST2toST2Type& Ct, const ST2toST2Type2& C, const TensorType& F) { - const auto C0 = F[0] * F[0]; - const auto C1 = F[1] * F[1]; - const auto C2 = F[2] * F[2]; - Ct(0, 0) = C0 * C0 * C(0, 0); - Ct(0, 1) = C0 * C1 * C(0, 1); - Ct(0, 2) = C0 * C2 * C(0, 2); - Ct(1, 0) = C1 * C0 * C(1, 0); - Ct(1, 1) = C1 * C1 * C(1, 1); - Ct(1, 2) = C1 * C2 * C(1, 2); - Ct(2, 0) = C2 * C0 * C(2, 0); - Ct(2, 1) = C2 * C1 * C(2, 1); - Ct(2, 2) = C2 * C2 * C(2, 2); - } - - template - typename std::enable_if() && - implementsST2toST2Concept() && - getSpaceDimension() == 2u && - getSpaceDimension() == 2u && - getSpaceDimension() == 2u, - void>::type - push_forward(ST2toST2Type& Ct, const ST2toST2Type2& C, const TensorType& F) { - using NumType = numeric_type; - constexpr auto cste = Cste::sqrt2; - constexpr auto icste = Cste::isqrt2; - Ct(0, 0) = F[0] * F[0] * F[0] * F[0] * C(0, 0) + - F[0] * F[0] * F[0] * F[3] * C(0, 3) * icste + - F[0] * F[0] * F[3] * F[0] * C(0, 3) * icste + - F[0] * F[0] * F[3] * F[3] * C(0, 1) + - F[0] * F[3] * F[0] * F[0] * C(3, 0) * icste + - F[0] * F[3] * F[0] * F[3] * C(3, 3) / 2 + - F[0] * F[3] * F[3] * F[0] * C(3, 3) / 2 + - F[0] * F[3] * F[3] * F[3] * C(3, 1) * icste + - F[3] * F[0] * F[0] * F[0] * C(3, 0) * icste + - F[3] * F[0] * F[0] * F[3] * C(3, 3) / 2 + - F[3] * F[0] * F[3] * F[0] * C(3, 3) / 2 + - F[3] * F[0] * F[3] * F[3] * C(3, 1) * icste + - F[3] * F[3] * F[0] * F[0] * C(1, 0) + - F[3] * F[3] * F[0] * F[3] * C(1, 3) * icste + - F[3] * F[3] * F[3] * F[0] * C(1, 3) * icste + - F[3] * F[3] * F[3] * F[3] * C(1, 1); - Ct(0, 3) = F[0] * F[0] * F[0] * F[4] * C(0, 0) * cste + - F[0] * F[0] * F[0] * F[1] * C(0, 3) + - F[0] * F[0] * F[3] * F[4] * C(0, 3) + - F[0] * F[0] * F[3] * F[1] * C(0, 1) * cste + - F[0] * F[3] * F[0] * F[4] * C(3, 0) + - F[0] * F[3] * F[0] * F[1] * C(3, 3) * icste + - F[0] * F[3] * F[3] * F[4] * C(3, 3) * icste + - F[0] * F[3] * F[3] * F[1] * C(3, 1) + - F[3] * F[0] * F[0] * F[4] * C(3, 0) + - F[3] * F[0] * F[0] * F[1] * C(3, 3) * icste + - F[3] * F[0] * F[3] * F[4] * C(3, 3) * icste + - F[3] * F[0] * F[3] * F[1] * C(3, 1) + - F[3] * F[3] * F[0] * F[4] * C(1, 0) * cste + - F[3] * F[3] * F[0] * F[1] * C(1, 3) + - F[3] * F[3] * F[3] * F[4] * C(1, 3) + - F[3] * F[3] * F[3] * F[1] * C(1, 1) * cste; - Ct(0, 1) = F[0] * F[0] * F[4] * F[4] * C(0, 0) + - F[0] * F[0] * F[4] * F[1] * C(0, 3) * icste + - F[0] * F[0] * F[1] * F[4] * C(0, 3) * icste + - F[0] * F[0] * F[1] * F[1] * C(0, 1) + - F[0] * F[3] * F[4] * F[4] * C(3, 0) * icste + - F[0] * F[3] * F[4] * F[1] * C(3, 3) / 2 + - F[0] * F[3] * F[1] * F[4] * C(3, 3) / 2 + - F[0] * F[3] * F[1] * F[1] * C(3, 1) * icste + - F[3] * F[0] * F[4] * F[4] * C(3, 0) * icste + - F[3] * F[0] * F[4] * F[1] * C(3, 3) / 2 + - F[3] * F[0] * F[1] * F[4] * C(3, 3) / 2 + - F[3] * F[0] * F[1] * F[1] * C(3, 1) * icste + - F[3] * F[3] * F[4] * F[4] * C(1, 0) + - F[3] * F[3] * F[4] * F[1] * C(1, 3) * icste + - F[3] * F[3] * F[1] * F[4] * C(1, 3) * icste + - F[3] * F[3] * F[1] * F[1] * C(1, 1); - Ct(0, 2) = F[0] * F[0] * F[2] * F[2] * C(0, 2) + - F[0] * F[3] * F[2] * F[2] * C(3, 2) * icste + - F[3] * F[0] * F[2] * F[2] * C(3, 2) * icste + - F[3] * F[3] * F[2] * F[2] * C(1, 2); - Ct(3, 0) = F[0] * F[4] * F[0] * F[0] * C(0, 0) * cste + - F[0] * F[4] * F[0] * F[3] * C(0, 3) + - F[0] * F[4] * F[3] * F[0] * C(0, 3) + - F[0] * F[4] * F[3] * F[3] * C(0, 1) * cste + - F[0] * F[1] * F[0] * F[0] * C(3, 0) + - F[0] * F[1] * F[0] * F[3] * C(3, 3) * icste + - F[0] * F[1] * F[3] * F[0] * C(3, 3) * icste + - F[0] * F[1] * F[3] * F[3] * C(3, 1) + - F[3] * F[4] * F[0] * F[0] * C(3, 0) + - F[3] * F[4] * F[0] * F[3] * C(3, 3) * icste + - F[3] * F[4] * F[3] * F[0] * C(3, 3) * icste + - F[3] * F[4] * F[3] * F[3] * C(3, 1) + - F[3] * F[1] * F[0] * F[0] * C(1, 0) * cste + - F[3] * F[1] * F[0] * F[3] * C(1, 3) + - F[3] * F[1] * F[3] * F[0] * C(1, 3) + - F[3] * F[1] * F[3] * F[3] * C(1, 1) * cste; - Ct(3, 3) = F[0] * F[4] * F[0] * F[4] * C(0, 0) * 2 + - F[0] * F[4] * F[0] * F[1] * C(0, 3) * cste + - F[0] * F[4] * F[3] * F[4] * C(0, 3) * cste + - F[0] * F[4] * F[3] * F[1] * C(0, 1) * 2 + - F[0] * F[1] * F[0] * F[4] * C(3, 0) * cste + - F[0] * F[1] * F[0] * F[1] * C(3, 3) + - F[0] * F[1] * F[3] * F[4] * C(3, 3) + - F[0] * F[1] * F[3] * F[1] * C(3, 1) * cste + - F[3] * F[4] * F[0] * F[4] * C(3, 0) * cste + - F[3] * F[4] * F[0] * F[1] * C(3, 3) + - F[3] * F[4] * F[3] * F[4] * C(3, 3) + - F[3] * F[4] * F[3] * F[1] * C(3, 1) * cste + - F[3] * F[1] * F[0] * F[4] * C(1, 0) * 2 + - F[3] * F[1] * F[0] * F[1] * C(1, 3) * cste + - F[3] * F[1] * F[3] * F[4] * C(1, 3) * cste + - F[3] * F[1] * F[3] * F[1] * C(1, 1) * 2; - Ct(3, 1) = F[0] * F[4] * F[4] * F[4] * C(0, 0) * cste + - F[0] * F[4] * F[4] * F[1] * C(0, 3) + - F[0] * F[4] * F[1] * F[4] * C(0, 3) + - F[0] * F[4] * F[1] * F[1] * C(0, 1) * cste + - F[0] * F[1] * F[4] * F[4] * C(3, 0) + - F[0] * F[1] * F[4] * F[1] * C(3, 3) * icste + - F[0] * F[1] * F[1] * F[4] * C(3, 3) * icste + - F[0] * F[1] * F[1] * F[1] * C(3, 1) + - F[3] * F[4] * F[4] * F[4] * C(3, 0) + - F[3] * F[4] * F[4] * F[1] * C(3, 3) * icste + - F[3] * F[4] * F[1] * F[4] * C(3, 3) * icste + - F[3] * F[4] * F[1] * F[1] * C(3, 1) + - F[3] * F[1] * F[4] * F[4] * C(1, 0) * cste + - F[3] * F[1] * F[4] * F[1] * C(1, 3) + - F[3] * F[1] * F[1] * F[4] * C(1, 3) + - F[3] * F[1] * F[1] * F[1] * C(1, 1) * cste; - Ct(3, 2) = F[0] * F[4] * F[2] * F[2] * C(0, 2) * cste + - F[0] * F[1] * F[2] * F[2] * C(3, 2) + - F[3] * F[4] * F[2] * F[2] * C(3, 2) + - F[3] * F[1] * F[2] * F[2] * C(1, 2) * cste; - Ct(1, 0) = F[4] * F[4] * F[0] * F[0] * C(0, 0) + - F[4] * F[4] * F[0] * F[3] * C(0, 3) * icste + - F[4] * F[4] * F[3] * F[0] * C(0, 3) * icste + - F[4] * F[4] * F[3] * F[3] * C(0, 1) + - F[4] * F[1] * F[0] * F[0] * C(3, 0) * icste + - F[4] * F[1] * F[0] * F[3] * C(3, 3) / 2 + - F[4] * F[1] * F[3] * F[0] * C(3, 3) / 2 + - F[4] * F[1] * F[3] * F[3] * C(3, 1) * icste + - F[1] * F[4] * F[0] * F[0] * C(3, 0) * icste + - F[1] * F[4] * F[0] * F[3] * C(3, 3) / 2 + - F[1] * F[4] * F[3] * F[0] * C(3, 3) / 2 + - F[1] * F[4] * F[3] * F[3] * C(3, 1) * icste + - F[1] * F[1] * F[0] * F[0] * C(1, 0) + - F[1] * F[1] * F[0] * F[3] * C(1, 3) * icste + - F[1] * F[1] * F[3] * F[0] * C(1, 3) * icste + - F[1] * F[1] * F[3] * F[3] * C(1, 1); - Ct(1, 3) = F[4] * F[4] * F[0] * F[4] * C(0, 0) * cste + - F[4] * F[4] * F[0] * F[1] * C(0, 3) + - F[4] * F[4] * F[3] * F[4] * C(0, 3) + - F[4] * F[4] * F[3] * F[1] * C(0, 1) * cste + - F[4] * F[1] * F[0] * F[4] * C(3, 0) + - F[4] * F[1] * F[0] * F[1] * C(3, 3) * icste + - F[4] * F[1] * F[3] * F[4] * C(3, 3) * icste + - F[4] * F[1] * F[3] * F[1] * C(3, 1) + - F[1] * F[4] * F[0] * F[4] * C(3, 0) + - F[1] * F[4] * F[0] * F[1] * C(3, 3) * icste + - F[1] * F[4] * F[3] * F[4] * C(3, 3) * icste + - F[1] * F[4] * F[3] * F[1] * C(3, 1) + - F[1] * F[1] * F[0] * F[4] * C(1, 0) * cste + - F[1] * F[1] * F[0] * F[1] * C(1, 3) + - F[1] * F[1] * F[3] * F[4] * C(1, 3) + - F[1] * F[1] * F[3] * F[1] * C(1, 1) * cste; - Ct(1, 1) = F[4] * F[4] * F[4] * F[4] * C(0, 0) + - F[4] * F[4] * F[4] * F[1] * C(0, 3) * icste + - F[4] * F[4] * F[1] * F[4] * C(0, 3) * icste + - F[4] * F[4] * F[1] * F[1] * C(0, 1) + - F[4] * F[1] * F[4] * F[4] * C(3, 0) * icste + - F[4] * F[1] * F[4] * F[1] * C(3, 3) / 2 + - F[4] * F[1] * F[1] * F[4] * C(3, 3) / 2 + - F[4] * F[1] * F[1] * F[1] * C(3, 1) * icste + - F[1] * F[4] * F[4] * F[4] * C(3, 0) * icste + - F[1] * F[4] * F[4] * F[1] * C(3, 3) / 2 + - F[1] * F[4] * F[1] * F[4] * C(3, 3) / 2 + - F[1] * F[4] * F[1] * F[1] * C(3, 1) * icste + - F[1] * F[1] * F[4] * F[4] * C(1, 0) + - F[1] * F[1] * F[4] * F[1] * C(1, 3) * icste + - F[1] * F[1] * F[1] * F[4] * C(1, 3) * icste + - F[1] * F[1] * F[1] * F[1] * C(1, 1); - Ct(1, 2) = F[4] * F[4] * F[2] * F[2] * C(0, 2) + - F[4] * F[1] * F[2] * F[2] * C(3, 2) * icste + - F[1] * F[4] * F[2] * F[2] * C(3, 2) * icste + - F[1] * F[1] * F[2] * F[2] * C(1, 2); - Ct(2, 0) = F[2] * F[2] * F[0] * F[0] * C(2, 0) + - F[2] * F[2] * F[0] * F[3] * C(2, 3) * icste + - F[2] * F[2] * F[3] * F[0] * C(2, 3) * icste + - F[2] * F[2] * F[3] * F[3] * C(2, 1); - Ct(2, 3) = F[2] * F[2] * F[0] * F[4] * C(2, 0) * cste + - F[2] * F[2] * F[0] * F[1] * C(2, 3) + - F[2] * F[2] * F[3] * F[4] * C(2, 3) + - F[2] * F[2] * F[3] * F[1] * C(2, 1) * cste; - Ct(2, 1) = F[2] * F[2] * F[4] * F[4] * C(2, 0) + - F[2] * F[2] * F[4] * F[1] * C(2, 3) * icste + - F[2] * F[2] * F[1] * F[4] * C(2, 3) * icste + - F[2] * F[2] * F[1] * F[1] * C(2, 1); - Ct(2, 2) = F[2] * F[2] * F[2] * F[2] * C(2, 2); - } - - template - typename std::enable_if() && - implementsST2toST2Concept() && - getSpaceDimension() == 3u && - getSpaceDimension() == 3u && - getSpaceDimension() == 3u, - void>::type - push_forward(ST2toST2Type& Ct, const ST2toST2Type2& C, const TensorType& F) { - using NumType = numeric_type; - using size_type = unsigned short; - auto row_index = [](size_type i, size_type j) -> size_type { - // i,j are valid for the space dimension considered - if ((i == j) && (i < 3)) { - return i; - } else if ((i == 0) && (j == 1)) { - return 3; - } else if ((i == 1) && (j == 0)) { - return 3; - } else if ((i == 0) && (j == 2)) { - return 4; - } else if ((i == 2) && (j == 0)) { - return 4; - } else if ((i == 1) && (j == 2)) { - return 5; - } - return 5; - }; - auto row_index2 = [](size_type i, size_type j) -> size_type { - // i,j are valid for the space dimension considered - if ((i == j) && (i < 3)) { - return i; - } else if ((i == 0) && (j == 1)) { - return 3; - } else if ((i == 1) && (j == 0)) { - return 4; - } else if ((i == 0) && (j == 2)) { - return 5; - } else if ((i == 2) && (j == 0)) { - return 6; - } else if ((i == 1) && (j == 2)) { - return 7; - } - return 8; - }; - auto set = [&Ct](const size_type i, const size_type j, const NumType v) { + TFEL_HOST_DEVICE constexpr void push_forward(ST2toST2Type& Ct, + const ST2toST2Type2& C, + const TensorType& F) noexcept // + requires(getSpaceDimension() == + getSpaceDimension() && + getSpaceDimension() == + getSpaceDimension()) { + constexpr auto N = getSpaceDimension(); + if constexpr (N == 1) { + const auto C0 = F[0] * F[0]; + const auto C1 = F[1] * F[1]; + const auto C2 = F[2] * F[2]; + Ct(0, 0) = C0 * C0 * C(0, 0); + Ct(0, 1) = C0 * C1 * C(0, 1); + Ct(0, 2) = C0 * C2 * C(0, 2); + Ct(1, 0) = C1 * C0 * C(1, 0); + Ct(1, 1) = C1 * C1 * C(1, 1); + Ct(1, 2) = C1 * C2 * C(1, 2); + Ct(2, 0) = C2 * C0 * C(2, 0); + Ct(2, 1) = C2 * C1 * C(2, 1); + Ct(2, 2) = C2 * C2 * C(2, 2); + } else if constexpr (N == 2) { + using NumType = numeric_type; constexpr auto cste = Cste::sqrt2; - if (((i > 2) && (j <= 2)) || ((j > 2) && (i <= 2))) { - Ct(i, j) = v * cste; - } else if ((i > 2) && (j > 2)) { - Ct(i, j) = v * 2; - } else { - Ct(i, j) = v; - } - }; - auto get = [&C](const size_type i, const size_type j) { constexpr auto icste = Cste::isqrt2; - if (((i > 2) && (j <= 2)) || ((j > 2) && (i <= 2))) { - return C(i, j) * icste; - } else if ((i > 2) && (j > 2)) { - return C(i, j) / 2; - } - return C(i, j); - }; - for (size_type i = 0; i != 3; ++i) { - for (size_type j = 0; j != 3; ++j) { - if (((i == 1) && (j == 0)) || ((i == 2) && (j == 0)) || - ((i == 2) && (j == 1))) { - continue; + Ct(0, 0) = F[0] * F[0] * F[0] * F[0] * C(0, 0) + + F[0] * F[0] * F[0] * F[3] * C(0, 3) * icste + + F[0] * F[0] * F[3] * F[0] * C(0, 3) * icste + + F[0] * F[0] * F[3] * F[3] * C(0, 1) + + F[0] * F[3] * F[0] * F[0] * C(3, 0) * icste + + F[0] * F[3] * F[0] * F[3] * C(3, 3) / 2 + + F[0] * F[3] * F[3] * F[0] * C(3, 3) / 2 + + F[0] * F[3] * F[3] * F[3] * C(3, 1) * icste + + F[3] * F[0] * F[0] * F[0] * C(3, 0) * icste + + F[3] * F[0] * F[0] * F[3] * C(3, 3) / 2 + + F[3] * F[0] * F[3] * F[0] * C(3, 3) / 2 + + F[3] * F[0] * F[3] * F[3] * C(3, 1) * icste + + F[3] * F[3] * F[0] * F[0] * C(1, 0) + + F[3] * F[3] * F[0] * F[3] * C(1, 3) * icste + + F[3] * F[3] * F[3] * F[0] * C(1, 3) * icste + + F[3] * F[3] * F[3] * F[3] * C(1, 1); + Ct(0, 3) = F[0] * F[0] * F[0] * F[4] * C(0, 0) * cste + + F[0] * F[0] * F[0] * F[1] * C(0, 3) + + F[0] * F[0] * F[3] * F[4] * C(0, 3) + + F[0] * F[0] * F[3] * F[1] * C(0, 1) * cste + + F[0] * F[3] * F[0] * F[4] * C(3, 0) + + F[0] * F[3] * F[0] * F[1] * C(3, 3) * icste + + F[0] * F[3] * F[3] * F[4] * C(3, 3) * icste + + F[0] * F[3] * F[3] * F[1] * C(3, 1) + + F[3] * F[0] * F[0] * F[4] * C(3, 0) + + F[3] * F[0] * F[0] * F[1] * C(3, 3) * icste + + F[3] * F[0] * F[3] * F[4] * C(3, 3) * icste + + F[3] * F[0] * F[3] * F[1] * C(3, 1) + + F[3] * F[3] * F[0] * F[4] * C(1, 0) * cste + + F[3] * F[3] * F[0] * F[1] * C(1, 3) + + F[3] * F[3] * F[3] * F[4] * C(1, 3) + + F[3] * F[3] * F[3] * F[1] * C(1, 1) * cste; + Ct(0, 1) = F[0] * F[0] * F[4] * F[4] * C(0, 0) + + F[0] * F[0] * F[4] * F[1] * C(0, 3) * icste + + F[0] * F[0] * F[1] * F[4] * C(0, 3) * icste + + F[0] * F[0] * F[1] * F[1] * C(0, 1) + + F[0] * F[3] * F[4] * F[4] * C(3, 0) * icste + + F[0] * F[3] * F[4] * F[1] * C(3, 3) / 2 + + F[0] * F[3] * F[1] * F[4] * C(3, 3) / 2 + + F[0] * F[3] * F[1] * F[1] * C(3, 1) * icste + + F[3] * F[0] * F[4] * F[4] * C(3, 0) * icste + + F[3] * F[0] * F[4] * F[1] * C(3, 3) / 2 + + F[3] * F[0] * F[1] * F[4] * C(3, 3) / 2 + + F[3] * F[0] * F[1] * F[1] * C(3, 1) * icste + + F[3] * F[3] * F[4] * F[4] * C(1, 0) + + F[3] * F[3] * F[4] * F[1] * C(1, 3) * icste + + F[3] * F[3] * F[1] * F[4] * C(1, 3) * icste + + F[3] * F[3] * F[1] * F[1] * C(1, 1); + Ct(0, 2) = F[0] * F[0] * F[2] * F[2] * C(0, 2) + + F[0] * F[3] * F[2] * F[2] * C(3, 2) * icste + + F[3] * F[0] * F[2] * F[2] * C(3, 2) * icste + + F[3] * F[3] * F[2] * F[2] * C(1, 2); + Ct(3, 0) = F[0] * F[4] * F[0] * F[0] * C(0, 0) * cste + + F[0] * F[4] * F[0] * F[3] * C(0, 3) + + F[0] * F[4] * F[3] * F[0] * C(0, 3) + + F[0] * F[4] * F[3] * F[3] * C(0, 1) * cste + + F[0] * F[1] * F[0] * F[0] * C(3, 0) + + F[0] * F[1] * F[0] * F[3] * C(3, 3) * icste + + F[0] * F[1] * F[3] * F[0] * C(3, 3) * icste + + F[0] * F[1] * F[3] * F[3] * C(3, 1) + + F[3] * F[4] * F[0] * F[0] * C(3, 0) + + F[3] * F[4] * F[0] * F[3] * C(3, 3) * icste + + F[3] * F[4] * F[3] * F[0] * C(3, 3) * icste + + F[3] * F[4] * F[3] * F[3] * C(3, 1) + + F[3] * F[1] * F[0] * F[0] * C(1, 0) * cste + + F[3] * F[1] * F[0] * F[3] * C(1, 3) + + F[3] * F[1] * F[3] * F[0] * C(1, 3) + + F[3] * F[1] * F[3] * F[3] * C(1, 1) * cste; + Ct(3, 3) = F[0] * F[4] * F[0] * F[4] * C(0, 0) * 2 + + F[0] * F[4] * F[0] * F[1] * C(0, 3) * cste + + F[0] * F[4] * F[3] * F[4] * C(0, 3) * cste + + F[0] * F[4] * F[3] * F[1] * C(0, 1) * 2 + + F[0] * F[1] * F[0] * F[4] * C(3, 0) * cste + + F[0] * F[1] * F[0] * F[1] * C(3, 3) + + F[0] * F[1] * F[3] * F[4] * C(3, 3) + + F[0] * F[1] * F[3] * F[1] * C(3, 1) * cste + + F[3] * F[4] * F[0] * F[4] * C(3, 0) * cste + + F[3] * F[4] * F[0] * F[1] * C(3, 3) + + F[3] * F[4] * F[3] * F[4] * C(3, 3) + + F[3] * F[4] * F[3] * F[1] * C(3, 1) * cste + + F[3] * F[1] * F[0] * F[4] * C(1, 0) * 2 + + F[3] * F[1] * F[0] * F[1] * C(1, 3) * cste + + F[3] * F[1] * F[3] * F[4] * C(1, 3) * cste + + F[3] * F[1] * F[3] * F[1] * C(1, 1) * 2; + Ct(3, 1) = F[0] * F[4] * F[4] * F[4] * C(0, 0) * cste + + F[0] * F[4] * F[4] * F[1] * C(0, 3) + + F[0] * F[4] * F[1] * F[4] * C(0, 3) + + F[0] * F[4] * F[1] * F[1] * C(0, 1) * cste + + F[0] * F[1] * F[4] * F[4] * C(3, 0) + + F[0] * F[1] * F[4] * F[1] * C(3, 3) * icste + + F[0] * F[1] * F[1] * F[4] * C(3, 3) * icste + + F[0] * F[1] * F[1] * F[1] * C(3, 1) + + F[3] * F[4] * F[4] * F[4] * C(3, 0) + + F[3] * F[4] * F[4] * F[1] * C(3, 3) * icste + + F[3] * F[4] * F[1] * F[4] * C(3, 3) * icste + + F[3] * F[4] * F[1] * F[1] * C(3, 1) + + F[3] * F[1] * F[4] * F[4] * C(1, 0) * cste + + F[3] * F[1] * F[4] * F[1] * C(1, 3) + + F[3] * F[1] * F[1] * F[4] * C(1, 3) + + F[3] * F[1] * F[1] * F[1] * C(1, 1) * cste; + Ct(3, 2) = F[0] * F[4] * F[2] * F[2] * C(0, 2) * cste + + F[0] * F[1] * F[2] * F[2] * C(3, 2) + + F[3] * F[4] * F[2] * F[2] * C(3, 2) + + F[3] * F[1] * F[2] * F[2] * C(1, 2) * cste; + Ct(1, 0) = F[4] * F[4] * F[0] * F[0] * C(0, 0) + + F[4] * F[4] * F[0] * F[3] * C(0, 3) * icste + + F[4] * F[4] * F[3] * F[0] * C(0, 3) * icste + + F[4] * F[4] * F[3] * F[3] * C(0, 1) + + F[4] * F[1] * F[0] * F[0] * C(3, 0) * icste + + F[4] * F[1] * F[0] * F[3] * C(3, 3) / 2 + + F[4] * F[1] * F[3] * F[0] * C(3, 3) / 2 + + F[4] * F[1] * F[3] * F[3] * C(3, 1) * icste + + F[1] * F[4] * F[0] * F[0] * C(3, 0) * icste + + F[1] * F[4] * F[0] * F[3] * C(3, 3) / 2 + + F[1] * F[4] * F[3] * F[0] * C(3, 3) / 2 + + F[1] * F[4] * F[3] * F[3] * C(3, 1) * icste + + F[1] * F[1] * F[0] * F[0] * C(1, 0) + + F[1] * F[1] * F[0] * F[3] * C(1, 3) * icste + + F[1] * F[1] * F[3] * F[0] * C(1, 3) * icste + + F[1] * F[1] * F[3] * F[3] * C(1, 1); + Ct(1, 3) = F[4] * F[4] * F[0] * F[4] * C(0, 0) * cste + + F[4] * F[4] * F[0] * F[1] * C(0, 3) + + F[4] * F[4] * F[3] * F[4] * C(0, 3) + + F[4] * F[4] * F[3] * F[1] * C(0, 1) * cste + + F[4] * F[1] * F[0] * F[4] * C(3, 0) + + F[4] * F[1] * F[0] * F[1] * C(3, 3) * icste + + F[4] * F[1] * F[3] * F[4] * C(3, 3) * icste + + F[4] * F[1] * F[3] * F[1] * C(3, 1) + + F[1] * F[4] * F[0] * F[4] * C(3, 0) + + F[1] * F[4] * F[0] * F[1] * C(3, 3) * icste + + F[1] * F[4] * F[3] * F[4] * C(3, 3) * icste + + F[1] * F[4] * F[3] * F[1] * C(3, 1) + + F[1] * F[1] * F[0] * F[4] * C(1, 0) * cste + + F[1] * F[1] * F[0] * F[1] * C(1, 3) + + F[1] * F[1] * F[3] * F[4] * C(1, 3) + + F[1] * F[1] * F[3] * F[1] * C(1, 1) * cste; + Ct(1, 1) = F[4] * F[4] * F[4] * F[4] * C(0, 0) + + F[4] * F[4] * F[4] * F[1] * C(0, 3) * icste + + F[4] * F[4] * F[1] * F[4] * C(0, 3) * icste + + F[4] * F[4] * F[1] * F[1] * C(0, 1) + + F[4] * F[1] * F[4] * F[4] * C(3, 0) * icste + + F[4] * F[1] * F[4] * F[1] * C(3, 3) / 2 + + F[4] * F[1] * F[1] * F[4] * C(3, 3) / 2 + + F[4] * F[1] * F[1] * F[1] * C(3, 1) * icste + + F[1] * F[4] * F[4] * F[4] * C(3, 0) * icste + + F[1] * F[4] * F[4] * F[1] * C(3, 3) / 2 + + F[1] * F[4] * F[1] * F[4] * C(3, 3) / 2 + + F[1] * F[4] * F[1] * F[1] * C(3, 1) * icste + + F[1] * F[1] * F[4] * F[4] * C(1, 0) + + F[1] * F[1] * F[4] * F[1] * C(1, 3) * icste + + F[1] * F[1] * F[1] * F[4] * C(1, 3) * icste + + F[1] * F[1] * F[1] * F[1] * C(1, 1); + Ct(1, 2) = F[4] * F[4] * F[2] * F[2] * C(0, 2) + + F[4] * F[1] * F[2] * F[2] * C(3, 2) * icste + + F[1] * F[4] * F[2] * F[2] * C(3, 2) * icste + + F[1] * F[1] * F[2] * F[2] * C(1, 2); + Ct(2, 0) = F[2] * F[2] * F[0] * F[0] * C(2, 0) + + F[2] * F[2] * F[0] * F[3] * C(2, 3) * icste + + F[2] * F[2] * F[3] * F[0] * C(2, 3) * icste + + F[2] * F[2] * F[3] * F[3] * C(2, 1); + Ct(2, 3) = F[2] * F[2] * F[0] * F[4] * C(2, 0) * cste + + F[2] * F[2] * F[0] * F[1] * C(2, 3) + + F[2] * F[2] * F[3] * F[4] * C(2, 3) + + F[2] * F[2] * F[3] * F[1] * C(2, 1) * cste; + Ct(2, 1) = F[2] * F[2] * F[4] * F[4] * C(2, 0) + + F[2] * F[2] * F[4] * F[1] * C(2, 3) * icste + + F[2] * F[2] * F[1] * F[4] * C(2, 3) * icste + + F[2] * F[2] * F[1] * F[1] * C(2, 1); + Ct(2, 2) = F[2] * F[2] * F[2] * F[2] * C(2, 2); + } else { + using NumType = numeric_type; + using size_type = unsigned short; + auto row_index = [](size_type i, size_type j) -> size_type { + // i,j are valid for the space dimension considered + if ((i == j) && (i < 3)) { + return i; + } else if ((i == 0) && (j == 1)) { + return 3; + } else if ((i == 1) && (j == 0)) { + return 3; + } else if ((i == 0) && (j == 2)) { + return 4; + } else if ((i == 2) && (j == 0)) { + return 4; + } else if ((i == 1) && (j == 2)) { + return 5; + } + return 5; + }; + auto row_index2 = [](size_type i, size_type j) -> size_type { + // i,j are valid for the space dimension considered + if ((i == j) && (i < 3)) { + return i; + } else if ((i == 0) && (j == 1)) { + return 3; + } else if ((i == 1) && (j == 0)) { + return 4; + } else if ((i == 0) && (j == 2)) { + return 5; + } else if ((i == 2) && (j == 0)) { + return 6; + } else if ((i == 1) && (j == 2)) { + return 7; } - const auto rij = row_index(i, j); - for (size_type k = 0; k != 3; ++k) { - for (size_type l = 0; l != 3; ++l) { - const auto rkl = row_index(k, l); - auto Cv = NumType(0); - for (size_type m = 0; m != 3; ++m) { - for (size_type n = 0; n != 3; ++n) { - for (size_type p = 0; p != 3; ++p) { - for (size_type q = 0; q != 3; ++q) { - const auto rim = row_index2(i, m); - const auto rjn = row_index2(j, n); - const auto rkp = row_index2(k, p); - const auto rlq = row_index2(l, q); - const auto rmn = row_index(m, n); - const auto rpq = row_index(p, q); - Cv += F[rim] * F[rjn] * F[rkp] * F[rlq] * get(rmn, rpq); + return 8; + }; + auto set = [&Ct](const size_type i, const size_type j, const NumType v) { + constexpr auto cste = Cste::sqrt2; + if (((i > 2) && (j <= 2)) || ((j > 2) && (i <= 2))) { + Ct(i, j) = v * cste; + } else if ((i > 2) && (j > 2)) { + Ct(i, j) = v * 2; + } else { + Ct(i, j) = v; + } + }; + auto get = [&C](const size_type i, const size_type j) { + constexpr auto icste = Cste::isqrt2; + if (((i > 2) && (j <= 2)) || ((j > 2) && (i <= 2))) { + return C(i, j) * icste; + } else if ((i > 2) && (j > 2)) { + return C(i, j) / 2; + } + return C(i, j); + }; + for (size_type i = 0; i != 3; ++i) { + for (size_type j = 0; j != 3; ++j) { + if (((i == 1) && (j == 0)) || ((i == 2) && (j == 0)) || + ((i == 2) && (j == 1))) { + continue; + } + const auto rij = row_index(i, j); + for (size_type k = 0; k != 3; ++k) { + for (size_type l = 0; l != 3; ++l) { + const auto rkl = row_index(k, l); + auto Cv = NumType(0); + for (size_type m = 0; m != 3; ++m) { + for (size_type n = 0; n != 3; ++n) { + for (size_type p = 0; p != 3; ++p) { + for (size_type q = 0; q != 3; ++q) { + const auto rim = row_index2(i, m); + const auto rjn = row_index2(j, n); + const auto rkp = row_index2(k, p); + const auto rlq = row_index2(l, q); + const auto rmn = row_index(m, n); + const auto rpq = row_index(p, q); + Cv += F[rim] * F[rjn] * F[rkp] * F[rlq] * get(rmn, rpq); + } } } } + set(rij, rkl, Cv); } - set(rij, rkl, Cv); } } } } } - } // end of namespace tfel::math #endif /* LIB_TFEL_MATH_ST2TOST2_PUSH_FORWARDIXX */ diff --git a/include/TFEL/Math/ST2toST2/ST2toST2ST2toST2ProductExpr.hxx b/include/TFEL/Math/ST2toST2/ST2toST2ST2toST2ProductExpr.hxx index 9ae9c4a41..9b3aa04a0 100644 --- a/include/TFEL/Math/ST2toST2/ST2toST2ST2toST2ProductExpr.hxx +++ b/include/TFEL/Math/ST2toST2/ST2toST2ST2toST2ProductExpr.hxx @@ -32,7 +32,7 @@ namespace tfel::math { template struct TFEL_VISIBILITY_LOCAL Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder< StensorDimeToSize()>::value * @@ -48,7 +48,7 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template + template TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a, const ST2toST2Type2& b) : array_holder()>::value, numeric_type>() { - static_assert(implementsST2toST2Concept()); - static_assert(implementsST2toST2Concept()); static_assert(getSpaceDimension() == 1u); static_assert(getSpaceDimension() == 1u); this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0); @@ -93,7 +91,7 @@ namespace tfel::math { template struct TFEL_VISIBILITY_LOCAL Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder< StensorDimeToSize()>::value * @@ -109,7 +107,7 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template + template TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a, const ST2toST2Type2& b) : array_holder()>::value, numeric_type>() { - static_assert(implementsST2toST2Concept()); - static_assert(implementsST2toST2Concept()); static_assert(getSpaceDimension() == 2u); static_assert(getSpaceDimension() == 2u); this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) + @@ -177,7 +173,7 @@ namespace tfel::math { template struct TFEL_VISIBILITY_LOCAL Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder< StensorDimeToSize()>::value * @@ -193,7 +189,7 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template + template TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a, const ST2toST2Type2& b) : array_holder()>::value, numeric_type>() { - static_assert(implementsST2toST2Concept()); - static_assert(implementsST2toST2Concept()); static_assert(getSpaceDimension() == 3u); static_assert(getSpaceDimension() == 3u); this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) + diff --git a/include/TFEL/Math/ST2toST2/ST2toST2StensorProductExpr.hxx b/include/TFEL/Math/ST2toST2/ST2toST2StensorProductExpr.hxx index 7e69d41de..af3603179 100644 --- a/include/TFEL/Math/ST2toST2/ST2toST2StensorProductExpr.hxx +++ b/include/TFEL/Math/ST2toST2/ST2toST2StensorProductExpr.hxx @@ -47,9 +47,9 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const ST2toST2Type& a, const StensorType& b) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a, + const StensorType& b) noexcept { static_assert(getSpaceDimension() == 1u); static_assert(getSpaceDimension() == 1u); this->v[0] = a(0, 0) * b(0) + a(0, 1) * b(1) + a(0, 2) * b(2); @@ -60,22 +60,20 @@ namespace tfel::math { * \brief access operator * \param[in] i : index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i) const noexcept { return this->v[i]; } // end of operator() /*! * \brief access operator * \param[in] i : index */ - TFEL_MATH_INLINE const value_type& operator[]( - const unsigned short i) const { + TFEL_HOST_DEVICE constexpr const value_type& operator[]( + const unsigned short i) const noexcept { return this->v[i]; } // end of operator[] - /*! - * \return the runtime properties - */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + //! \return the runtime properties + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr @@ -99,9 +97,9 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const ST2toST2Type& a, const StensorType& b) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a, + const StensorType& b) noexcept { static_assert(getSpaceDimension() == 2u); static_assert(getSpaceDimension() == 2u); this->v[0] = @@ -117,22 +115,21 @@ namespace tfel::math { * \brief access operator * \param[in] i : index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i) const noexcept { return this->v[i]; } // end of operator() /*! * \brief access operator * \param[in] i : index */ - TFEL_MATH_INLINE const value_type& operator[]( - const unsigned short i) const { + TFEL_HOST_DEVICE constexpr const value_type& operator[]( + const unsigned short i) const noexcept { return this->v[i]; } // end of operator[] - /*! - * \return the runtime properties - */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + //! \return the runtime properties + TFEL_HOST_DEVICE constexpr RunTimeProperties getRunTimeProperties() const + noexcept { return RunTimeProperties(); } }; // end of struct Expr @@ -156,9 +153,9 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const ST2toST2Type& a, const StensorType& b) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a, + const StensorType& b) noexcept { static_assert(getSpaceDimension() == 3u); static_assert(getSpaceDimension() == 3u); this->v[0] = a(0, 0) * b(0) + a(0, 1) * b(1) + a(0, 2) * b(2) + @@ -178,22 +175,20 @@ namespace tfel::math { * \brief access operator * \param[in] i : index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i) const noexcept { return this->v[i]; } // end of operator() /*! * \brief access operator * \param[in] i : index */ - TFEL_MATH_INLINE const value_type& operator[]( - const unsigned short i) const { + TFEL_HOST_DEVICE constexpr const value_type& operator[]( + const unsigned short i) const noexcept { return this->v[i]; } // end of operator[] - /*! - * \return the runtime properties - */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + //! \return the runtime properties + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr diff --git a/include/TFEL/Math/ST2toST2/ST2toST2TransposeExpr.hxx b/include/TFEL/Math/ST2toST2/ST2toST2TransposeExpr.hxx index b20458979..4c6834426 100644 --- a/include/TFEL/Math/ST2toST2/ST2toST2TransposeExpr.hxx +++ b/include/TFEL/Math/ST2toST2/ST2toST2TransposeExpr.hxx @@ -30,17 +30,16 @@ namespace tfel::math { * \tparam A: type of the object to be transposed * \pre A must match the `ST2toST2Concept` concept */ - template + template struct TFEL_VISIBILITY_LOCAL ST2toST2TransposeExpr : public ExprBase { - static_assert(implementsST2toST2Concept()); //! a simple alias using RunTimeProperties = EmptyRunTimeProperties; //! a simple alias using IndexType = unsigned short; //! a simple alias using NumType = numeric_type; - - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + //! \return the runtime properties + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return EmptyRunTimeProperties(); } @@ -52,11 +51,20 @@ namespace tfel::math { typedef const NumType& const_reference; typedef IndexType size_type; typedef ptrdiff_t difference_type; - - TFEL_MATH_INLINE ST2toST2TransposeExpr(A l) : a(l) {} - - TFEL_MATH_INLINE NumType operator()(const IndexType i, - const IndexType j) const { + /*! + * \brief access operator + * \param[in] i: row number + * \param[in] j: colum number + */ + TFEL_HOST_DEVICE constexpr ST2toST2TransposeExpr(A l) noexcept : a(l) {} + /*! + * \brief access operator + * \param[in] i: row number + * \param[in] j: colum number + */ + TFEL_HOST_DEVICE constexpr NumType operator()(const IndexType i, + const IndexType j) const + noexcept { return this->a(j, i); } // end of operator() //! storage for the object diff --git a/include/TFEL/Math/ST2toST2/StensorST2toST2ProductExpr.hxx b/include/TFEL/Math/ST2toST2/StensorST2toST2ProductExpr.hxx index 58500eacb..4083484b6 100644 --- a/include/TFEL/Math/ST2toST2/StensorST2toST2ProductExpr.hxx +++ b/include/TFEL/Math/ST2toST2/StensorST2toST2ProductExpr.hxx @@ -27,9 +27,7 @@ namespace tfel::math { struct StensorST2toST2ProductExpr { }; // end of struct StensorST2toST2ProductExpr - /*! - * Partial specialisation - */ + //!\brief partial specialisation in 1D template struct Expr> : public StensorConceptBase< @@ -47,9 +45,9 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const StensorType& a, const ST2toST2Type& b) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& a, + const ST2toST2Type& b) noexcept { static_assert(getSpaceDimension() == 1u); static_assert(getSpaceDimension() == 1u); this->v[0] = b(0, 0) * a(0) + b(1, 0) * a(1) + b(2, 0) * a(2); @@ -60,21 +58,25 @@ namespace tfel::math { * \brief access operator * \param[in] i : index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i) const { - return this->v[i]; - } // end of operator() + TFEL_HOST_DEVICE constexpr const value_type& operator[]( + const unsigned short i) const noexcept { + return this->operator()(i); + } /*! - * \return the runtime properties + * \brief access operator + * \param[in] i : index */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i) const noexcept { + return this->v[i]; + } // end of operator() + //! \return the runtime properties + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr - /*! - * Partial specialisation - */ + //! \brief partial specialisation in 2D template struct Expr> : public StensorConceptBase< @@ -91,9 +93,9 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const StensorType& a, const ST2toST2Type& b) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& a, + const ST2toST2Type& b) noexcept { static_assert(getSpaceDimension() == 2u); static_assert(getSpaceDimension() == 2u); this->v[0] = @@ -109,21 +111,25 @@ namespace tfel::math { * \brief access operator * \param[in] i : index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i) const { - return this->v[i]; - } // end of operator() + TFEL_HOST_DEVICE constexpr const value_type& operator[]( + const unsigned short i) const noexcept { + return this->operator()(i); + } /*! - * \return the runtime properties + * \brief access operator + * \param[in] i : index */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i) const noexcept { + return this->v[i]; + } // end of operator() + //! \return the runtime properties + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr - /*! - * Partial specialisation - */ + //! \brief partial specialisation in 3D template struct Expr> : public StensorConceptBase< @@ -140,9 +146,9 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const StensorType& a, const ST2toST2Type& b) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& a, + const ST2toST2Type& b) noexcept { static_assert(getSpaceDimension() == 3u); static_assert(getSpaceDimension() == 3u); this->v[0] = b(0, 0) * a(0) + b(1, 0) * a(1) + b(2, 0) * a(2) + @@ -162,14 +168,20 @@ namespace tfel::math { * \brief access operator * \param[in] i : index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i) const { - return this->v[i]; - } // end of operator() + TFEL_HOST_DEVICE constexpr const value_type& operator[]( + const unsigned short i) const noexcept { + return this->operator()(i); + } /*! - * \return the runtime properties + * \brief access operator + * \param[in] i : index */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i) const noexcept { + return this->v[i]; + } // end of operator() + //! \return the runtime properties + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr diff --git a/include/TFEL/Math/ST2toST2/StensorSquareDerivative.hxx b/include/TFEL/Math/ST2toST2/StensorSquareDerivative.hxx index aea5c8a48..5159628f7 100644 --- a/include/TFEL/Math/ST2toST2/StensorSquareDerivative.hxx +++ b/include/TFEL/Math/ST2toST2/StensorSquareDerivative.hxx @@ -19,19 +19,15 @@ namespace tfel::math { - /*! - * Empty structure allowing partial specialisation - */ + //! \brief empty structure allowing partial specialisation template struct StensorSquareDerivativeExpr { }; // end of struct StensorSquareDerivativeExpr - /*! - * Partial specialisation for 1D tensor - */ - template + //! \brief partial specialisation for 1D tensor + template struct Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder<9u, numeric_type> { static_assert(getSpaceDimension() == 1u); @@ -43,7 +39,7 @@ namespace tfel::math { * \param[in] B : second tensor of the product */ template - Expr(const StensorType& B) { + TFEL_HOST_DEVICE constexpr Expr(const StensorType& B) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(isAssignableTo, @@ -60,9 +56,9 @@ namespace tfel::math { * \param[in] B : second tensor of the product * \param[in] C : derivative of the first tensor */ - template - Expr(const StensorType& s, const ST2toST2Type& C) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& s, + const ST2toST2Type& C) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(getSpaceDimension() == @@ -87,8 +83,8 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - const value_type& operator()(const unsigned short i, - const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 3 + j]; } // end of operator() /*! @@ -96,17 +92,15 @@ namespace tfel::math { * In this case, the number of lines and columns * are deduced from the template parameter */ - RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr > - /*! - * Partial specialisation for 2D tensor - */ - template + //! \brief partial specialisation for 2D tensor + template struct Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder<16u, numeric_type> { static_assert(getSpaceDimension() == 2u); @@ -118,7 +112,7 @@ namespace tfel::math { * \param[in] B : second tensor of the product */ template - Expr(const StensorType& s) { + TFEL_HOST_DEVICE constexpr Expr(const StensorType& s) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(isAssignableTo, @@ -145,9 +139,9 @@ namespace tfel::math { * \param[in] B : second tensor of the product * \param[in] C : derivative of the first tensor */ - template - Expr(const StensorType& s, const ST2toST2Type& C) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& s, + const ST2toST2Type& C) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(getSpaceDimension() == @@ -179,8 +173,8 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - const value_type& operator()(const unsigned short i, - const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 4 + j]; } // end of operator() /*! @@ -188,17 +182,15 @@ namespace tfel::math { * In this case, the number of lines and columns * are deduced from the template parameter */ - RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr > - /*! - * Partial specialisation for 3D tensor - */ - template + //! \brief partial specialisation for 3D tensor + template struct Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder<36u, numeric_type> { static_assert(getSpaceDimension() == 3u); @@ -210,7 +202,7 @@ namespace tfel::math { * \param[in] B : second tensor of the product */ template - Expr(const StensorType& s) { + TFEL_HOST_DEVICE constexpr Expr(const StensorType& s) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(isAssignableTo, @@ -260,9 +252,9 @@ namespace tfel::math { * \param[in] B : second tensor of the product * \param[in] C : derivative of the first tensor */ - template - Expr(const StensorType& s, const ST2toST2Type& C) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& s, + const ST2toST2Type& C) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(getSpaceDimension() == @@ -371,8 +363,8 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - const value_type& operator()(const unsigned short i, - const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 6 + j]; } // end of operator() /*! @@ -380,7 +372,7 @@ namespace tfel::math { * In this case, the number of lines and columns * are deduced from the template parameter */ - RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr > diff --git a/include/TFEL/Math/ST2toST2/st2tost2.ixx b/include/TFEL/Math/ST2toST2/st2tost2.ixx index 7abb693b0..1bd809d32 100644 --- a/include/TFEL/Math/ST2toST2/st2tost2.ixx +++ b/include/TFEL/Math/ST2toST2/st2tost2.ixx @@ -26,33 +26,30 @@ #include "TFEL/Math/ST2toST2/BuildFromRotationMatrix.hxx" #include "TFEL/Math/ST2toST2/StensorSymmetricProductDerivative.hxx" #include "TFEL/Math/ST2toST2/SymmetricStensorProductDerivative.hxx" -#include "TFEL/Math/ST2toST2/ChangeBasis.hxx" namespace tfel::math { template template - std::enable_if_t() == N && - isAssignableTo, T>(), - Expr, StensorSquareDerivativeExpr>> - st2tost2::dsquare(const StensorType& s) { + TFEL_HOST_DEVICE constexpr auto + st2tost2::dsquare(const StensorType& s) noexcept requires( + getSpaceDimension() == N && + isAssignableTo, T>()) { return Expr, StensorSquareDerivativeExpr>(s); } template - template - std::enable_if_t< - implementsST2toST2Concept() && + template + TFEL_HOST_DEVICE constexpr auto st2tost2:: + dsquare(const StensorType& s, const ST2toST2Type& C) noexcept requires( getSpaceDimension() == N && getSpaceDimension() == N && isAssignableTo, numeric_type, OpMult>, - T>(), - Expr, StensorSquareDerivativeExpr>> - st2tost2::dsquare(const StensorType& s, const ST2toST2Type& C) { + T>()) { return Expr, StensorSquareDerivativeExpr>(s, C); - } + } // end of dsquare template template @@ -80,7 +77,7 @@ namespace tfel::math { } // end of st2tost2::fromRotationMatrix template - constexpr st2tost2 st2tost2::Id() { + constexpr st2tost2 st2tost2::Id() noexcept { constexpr auto c0 = T{0}; constexpr auto c1 = T{1}; static_assert((N == 1) || (N == 2) || (N == 3)); @@ -104,7 +101,7 @@ namespace tfel::math { } // end of st2tost2::Id template - constexpr st2tost2 st2tost2::IxI() { + constexpr st2tost2 st2tost2::IxI() noexcept { constexpr auto c1 = T{1}; static_assert((N == 1) || (N == 2) || (N == 3)); if constexpr (N == 1) { @@ -129,7 +126,7 @@ namespace tfel::math { } // end of st2tost2::Id template - constexpr st2tost2 st2tost2::K() { + constexpr st2tost2 st2tost2::K() noexcept { constexpr auto c2_3 = T{2} / T{3}; constexpr auto mc1_3 = -T{1} / T{3}; static_assert((N == 1) || (N == 2) || (N == 3)); @@ -157,7 +154,7 @@ namespace tfel::math { } // end of st2tost2::K template - constexpr st2tost2 st2tost2::M() { + constexpr st2tost2 st2tost2::M() noexcept { constexpr auto c1 = T{1}; constexpr auto mc1_2 = -T{1} / T{2}; static_assert((N == 1) || (N == 2) || (N == 3)); @@ -185,7 +182,7 @@ namespace tfel::math { } // end of st2tost2::M template - constexpr st2tost2 st2tost2::J() { + constexpr st2tost2 st2tost2::J() noexcept { constexpr auto c1_3 = T{1} / T{3}; static_assert((N == 1) || (N == 2) || (N == 3)); if constexpr (N == 1) { @@ -210,22 +207,24 @@ namespace tfel::math { } // end of st2tost2::J template - template - TFEL_MATH_INLINE2 void st2tost2::copy(const InputIterator src) { + TFEL_HOST_DEVICE constexpr void st2tost2::import( + const base_type* const src) noexcept { + tfel::fsalgo::transform< + StensorDimeToSize::value * + StensorDimeToSize::value>::exe(src, this->begin(), + [](const auto& v) { return T(v); }); + } + + template + TFEL_HOST_DEVICE constexpr void st2tost2::copy(const auto p) noexcept { tfel::fsalgo::copy::value * - StensorDimeToSize::value>::exe(src, *this); + StensorDimeToSize::value>::exe(p, *this); } - template - TFEL_MATH_INLINE2 std::enable_if_t< - implementsST2toST2Concept(), - st2tost2(), - BinaryOperationResult>, - numeric_type, - OpDiv>>> - invert(const ST2toST2Type& s) { - static constexpr unsigned short N = getSpaceDimension(); - static constexpr unsigned short StensorSize = StensorDimeToSize::value; + TFEL_HOST constexpr auto invert(const ST2toST2Concept auto& s) { + using ST2toST2Type = decltype(s); + constexpr auto N = getSpaceDimension(); + constexpr auto StensorSize = StensorDimeToSize::value; using NumType = numeric_type; using real = base_type; using iNumType = BinaryOperationResult; @@ -245,144 +244,142 @@ namespace tfel::math { return is; } // end of invert - template - std::enable_if_t< - implementsST2toST2Concept(), - st2tost2(), numeric_type>> + template + TFEL_HOST_DEVICE constexpr st2tost2(), + numeric_type> change_basis(const ST2toST2Type& s, - const rotation_matrix>& r) { - return st2tost2_internals::ChangeBasis< - getSpaceDimension()>::exe(s, r); - } + const rotation_matrix>& r) noexcept { + constexpr auto N = getSpaceDimension(); + if constexpr (N == 1) { + return s; + } else { + using st2tost2 = + tfel::math::st2tost2>; + const auto sr = st2tost2::fromRotationMatrix(r); + const auto sir = st2tost2::fromRotationMatrix(transpose(r)); + return sr * s * sir; + } + } // end of change_basie - template - std::enable_if_t() && - getSpaceDimension() == - getSpaceDimension(), - st2tost2(), - BinaryOperationResult, - numeric_type, - OpMult>>> - push_forward(const ST2toST2Type& C, const TensorType& F) { + template + TFEL_HOST_DEVICE constexpr auto + push_forward(const ST2toST2Type& C, const TensorType& F) noexcept requires( + getSpaceDimension() == getSpaceDimension()) { st2tost2(), BinaryOperationResult, numeric_type, OpMult>> r; push_forward(r, C, F); return r; - } + } // end of push_forward - template - std::enable_if_t< - implementsST2toST2Concept() && - getSpaceDimension() == getSpaceDimension(), - st2tost2(), - typename ComputeBinaryResult, - numeric_type, - OpMult>::Result>> - pull_back(const ST2toST2Type& C, const TensorType& F) { + template + TFEL_HOST constexpr auto pull_back( + const ST2toST2Type& C, + const TensorType& F) requires(getSpaceDimension() == + getSpaceDimension()) { const auto iF = invert(F); return push_forward(C, iF); - } + } // end of pull_back - template - std::enable_if_t< - isScalar>(), - st2tost2(), numeric_type>> - computeDeterminantSecondDerivative(const StensorType& s) { + TFEL_HOST_DEVICE constexpr auto computeDeterminantSecondDerivative( + const StensorConcept auto& s) noexcept { + using StensorType = decltype(s); using NumType = numeric_type; constexpr auto N = getSpaceDimension(); constexpr auto zero = NumType{0}; static_assert((N == 1) || (N == 2) || (N == 3)); + using Result = st2tost2>; if constexpr (N == 1) { - return {zero, s[2], s[1], s[2], zero, s[0], s[1], s[0], zero}; + return Result{zero, s[2], s[1], s[2], zero, s[0], s[1], s[0], zero}; } else if constexpr (N == 2) { - return {zero, s[2], s[1], zero, s[2], zero, s[0], zero, - s[1], s[0], zero, -s[3], zero, zero, -s[3], -s[2]}; + return Result{zero, s[2], s[1], zero, s[2], zero, s[0], zero, + s[1], s[0], zero, -s[3], zero, zero, -s[3], -s[2]}; } else { constexpr auto icste = Cste::isqrt2; - return {zero, s[2], s[1], zero, zero, -s[5], - s[2], zero, s[0], zero, -s[4], zero, - s[1], s[0], zero, -s[3], zero, zero, - zero, zero, -s[3], -s[2], s[5] * icste, s[4] * icste, - zero, -s[4], zero, s[5] * icste, -s[1], s[3] * icste, - -s[5], zero, zero, s[4] * icste, s[3] * icste, -s[0]}; + return Result{ + zero, s[2], s[1], zero, zero, -s[5], + s[2], zero, s[0], zero, -s[4], zero, + s[1], s[0], zero, -s[3], zero, zero, + zero, zero, -s[3], -s[2], s[5] * icste, s[4] * icste, + zero, -s[4], zero, s[5] * icste, -s[1], s[3] * icste, + -s[5], zero, zero, s[4] * icste, s[3] * icste, -s[0]}; } } // end of computeDeterminantSecondDerivative - template - std::enable_if_t< - isScalar>(), - st2tost2(), numeric_type>> - computeDeviatorDeterminantSecondDerivative(const StensorType& s) { + TFEL_HOST_DEVICE constexpr auto computeDeviatorDeterminantSecondDerivative( + const StensorConcept auto& s) noexcept { + using StensorType = decltype(s); constexpr auto N = getSpaceDimension(); static_assert((N == 1) || (N == 2) || (N == 3)); + using Result = + st2tost2(), numeric_type>; if constexpr (N == 1) { - return {-(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, - (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, - -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, - (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, - -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, - -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, - -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, - -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, - (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9}; + return Result{-(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, + (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, + -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, + (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, + -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, + -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, + -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, + -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, + (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9}; } else if constexpr (N == 2) { - return {-(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, - (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, - -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, - s[3] / 3, - (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, - -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, - -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, - s[3] / 3, - -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, - -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, - (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, - (-2 * s[3]) / 3, - s[3] / 3, - s[3] / 3, - (-2 * s[3]) / 3, - -(2 * s[2] - s[1] - s[0]) / 3}; + return Result{-(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, + (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, + -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, + s[3] / 3, + (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, + -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, + -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, + s[3] / 3, + -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, + -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, + (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, + (-2 * s[3]) / 3, + s[3] / 3, + s[3] / 3, + (-2 * s[3]) / 3, + -(2 * s[2] - s[1] - s[0]) / 3}; } else { using NumType = numeric_type; constexpr auto icste = Cste::isqrt2; - return {-(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, - (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, - -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, - s[3] / 3, - s[4] / 3, - (-2 * s[5]) / 3, - (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, - -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, - -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, - s[3] / 3, - (-2 * s[4]) / 3, - s[5] / 3, - -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, - -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, - (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, - (-2 * s[3]) / 3, - s[4] / 3, - s[5] / 3, - s[3] / 3, - s[3] / 3, - (-2 * s[3]) / 3, - -(2 * s[2] - s[1] - s[0]) / 3, - s[5] * icste, - s[4] * icste, - s[4] / 3, - (-2 * s[4]) / 3, - s[4] / 3, - s[5] * icste, - (s[2] - 2 * s[1] + s[0]) / 3, - s[3] * icste, - (-2 * s[5]) / 3, - s[5] / 3, - s[5] / 3, - s[4] * icste, - s[3] * icste, - (s[2] + s[1] - 2 * s[0]) / 3}; + return Result{-(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, + (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, + -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, + s[3] / 3, + s[4] / 3, + (-2 * s[5]) / 3, + (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, + -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, + -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, + s[3] / 3, + (-2 * s[4]) / 3, + s[5] / 3, + -(2 * s[2] - 4 * s[1] + 2 * s[0]) / 9, + -(2 * s[2] + 2 * s[1] - 4 * s[0]) / 9, + (4 * s[2] - 2 * s[1] - 2 * s[0]) / 9, + (-2 * s[3]) / 3, + s[4] / 3, + s[5] / 3, + s[3] / 3, + s[3] / 3, + (-2 * s[3]) / 3, + -(2 * s[2] - s[1] - s[0]) / 3, + s[5] * icste, + s[4] * icste, + s[4] / 3, + (-2 * s[4]) / 3, + s[4] / 3, + s[5] * icste, + (s[2] - 2 * s[1] + s[0]) / 3, + s[3] * icste, + (-2 * s[5]) / 3, + s[5] / 3, + s[5] / 3, + s[4] * icste, + s[3] * icste, + (s[2] + s[1] - 2 * s[0]) / 3}; } } // end of computeDeviatorDeterminantSecondDerivative diff --git a/include/TFEL/Math/ST2toT2/ST2toT2ST2toST2ProductExpr.hxx b/include/TFEL/Math/ST2toT2/ST2toT2ST2toST2ProductExpr.hxx index cc006e55b..ce90d04cf 100644 --- a/include/TFEL/Math/ST2toT2/ST2toT2ST2toST2ProductExpr.hxx +++ b/include/TFEL/Math/ST2toT2/ST2toT2ST2toST2ProductExpr.hxx @@ -50,10 +50,10 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const ST2toT2Type& a, const ST2toST2Type2& b) { + template + TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a, + const ST2toST2Type2& b) noexcept { static_assert(implementsST2toT2Concept()); - static_assert(implementsST2toST2Concept()); static_assert(getSpaceDimension() == 1u); static_assert(getSpaceDimension() == 1u); this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0); @@ -71,14 +71,14 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i, const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 3 + j]; } // end of operator() /*! * \return the runtime properties */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr @@ -105,10 +105,10 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const ST2toT2Type& a, const ST2toST2Type2& b) { + template + TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a, + const ST2toST2Type2& b) noexcept { static_assert(implementsST2toT2Concept()); - static_assert(implementsST2toST2Concept()); static_assert(getSpaceDimension() == 2u); static_assert(getSpaceDimension() == 2u); this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) + @@ -157,14 +157,14 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i, const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 4 + j]; } // end of operator() /*! * \return the runtime properties */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr @@ -191,10 +191,10 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const ST2toT2Type& a, const ST2toST2Type2& b) { + template + TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a, + const ST2toST2Type2& b) noexcept { static_assert(implementsST2toT2Concept()); - static_assert(implementsST2toST2Concept()); static_assert(getSpaceDimension() == 3u); static_assert(getSpaceDimension() == 3u); this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) + @@ -311,14 +311,14 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i, const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 6 + j]; } // end of operator() /*! * \return the runtime properties */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr diff --git a/include/TFEL/Math/ST2toT2/StensorProductLeftDerivativeExpr.hxx b/include/TFEL/Math/ST2toT2/StensorProductLeftDerivativeExpr.hxx index 839fc766d..24d0686a3 100644 --- a/include/TFEL/Math/ST2toT2/StensorProductLeftDerivativeExpr.hxx +++ b/include/TFEL/Math/ST2toT2/StensorProductLeftDerivativeExpr.hxx @@ -18,9 +18,7 @@ namespace tfel::math { - /*! - * Empty structure allowing partial specialisation - */ + //! \brief empty structure allowing partial specialisation template struct StensorProductLeftDerivativeExpr { }; // end of struct StensorProductLeftDerivativeExpr @@ -42,7 +40,7 @@ namespace tfel::math { * \param[in] B : second tensor of the product */ template - Expr(const StensorType& b) { + TFEL_HOST_DEVICE constexpr Expr(const StensorType& b) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(isAssignableTo, @@ -59,9 +57,9 @@ namespace tfel::math { * \param[in] B : second tensor of the product * \param[in] C : derivative of the first tensor */ - template - Expr(const StensorType& b, const ST2toST2Type& C) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& b, + const ST2toST2Type& C) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(getSpaceDimension() == @@ -86,8 +84,8 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - const value_type& operator()(const unsigned short i, - const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 3 + j]; } // end of operator() /*! @@ -95,7 +93,7 @@ namespace tfel::math { * In this case, the number of lines and columns * are deduced from the template parameter */ - RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct @@ -118,7 +116,7 @@ namespace tfel::math { * \param[in] B : second tensor of the product */ template - Expr(const StensorType& b) { + TFEL_HOST_DEVICE constexpr Expr(const StensorType& b) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(isAssignableTo, @@ -152,9 +150,9 @@ namespace tfel::math { * \param[in] B : second tensor of the product * \param[in] C : derivative of the first tensor */ - template - Expr(const StensorType& b, const ST2toST2Type& C) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& b, + const ST2toST2Type& C) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(getSpaceDimension() == @@ -193,8 +191,8 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - const value_type& operator()(const unsigned short i, - const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 4 + j]; } // end of operator() /*! @@ -202,11 +200,11 @@ namespace tfel::math { * In this case, the number of lines and columns * are deduced from the template parameter */ - RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct - // Expr > + // Expr > /*! * Partial specialisation for 3D tensor @@ -225,7 +223,7 @@ namespace tfel::math { * \param[in] B : second tensor of the product */ template - Expr(const StensorType& b) { + TFEL_HOST_DEVICE constexpr Expr(const StensorType& b) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(isAssignableTo, @@ -293,9 +291,9 @@ namespace tfel::math { * \param[in] B : second tensor of the product * \param[in] C : derivative of the first tensor */ - template - Expr(const StensorType& b, const ST2toST2Type& C) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& b, + const ST2toST2Type& C) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(getSpaceDimension() == @@ -421,8 +419,8 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - const value_type& operator()(const unsigned short i, - const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 6 + j]; } // end of operator() /*! @@ -430,11 +428,11 @@ namespace tfel::math { * In this case, the number of lines and columns * are deduced from the template parameter */ - RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct - // Expr > + // Expr > } // namespace tfel::math diff --git a/include/TFEL/Math/ST2toT2/StensorProductRightDerivativeExpr.hxx b/include/TFEL/Math/ST2toT2/StensorProductRightDerivativeExpr.hxx index 7ccfb3f57..b28847873 100644 --- a/include/TFEL/Math/ST2toT2/StensorProductRightDerivativeExpr.hxx +++ b/include/TFEL/Math/ST2toT2/StensorProductRightDerivativeExpr.hxx @@ -43,7 +43,7 @@ namespace tfel::math { * \param[in] a : first stensor of the product */ template - Expr(const StensorType& a) { + TFEL_HOST_DEVICE constexpr Expr(const StensorType& a) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(isAssignableTo, @@ -60,9 +60,9 @@ namespace tfel::math { * \param[in] A : first tensor of the product * \param[in] C : derivative of the second tensor */ - template - Expr(const StensorType& a, const ST2toST2Type& C) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& a, + const ST2toST2Type& C) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(getSpaceDimension() == @@ -87,8 +87,8 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - const value_type& operator()(const unsigned short i, - const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 3 + j]; } // end of operator() /*! @@ -96,7 +96,7 @@ namespace tfel::math { * In this case, the number of lines and columns * are deduced from the template parameter */ - RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct @@ -119,7 +119,7 @@ namespace tfel::math { * \param[in] a : first stensor of the product */ template - Expr(const StensorType& a) { + TFEL_HOST_DEVICE constexpr Expr(const StensorType& a) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(isAssignableTo, @@ -153,9 +153,9 @@ namespace tfel::math { * \param[in] A : first tensor of the product * \param[in] C : derivative of the second tensor */ - template - Expr(const StensorType& a, const ST2toST2Type& C) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& a, + const ST2toST2Type& C) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(getSpaceDimension() == @@ -194,8 +194,8 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - const value_type& operator()(const unsigned short i, - const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 4 + j]; } // end of operator() /*! @@ -203,11 +203,11 @@ namespace tfel::math { * In this case, the number of lines and columns * are deduced from the template parameter */ - RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct - // Expr > + // Expr > /*! * Partial specialisation for 3D stensor @@ -226,7 +226,7 @@ namespace tfel::math { * \param[in] a : first stensor of the product */ template - Expr(const StensorType& a) { + TFEL_HOST_DEVICE constexpr Expr(const StensorType& a) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(isAssignableTo, @@ -294,9 +294,9 @@ namespace tfel::math { * \param[in] a : first stensor of the product * \param[in] C : derivative of the second tensor */ - template - Expr(const StensorType& a, const ST2toST2Type& C) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const StensorType& a, + const ST2toST2Type& C) noexcept { static_assert(getSpaceDimension() == getSpaceDimension()); static_assert(getSpaceDimension() == @@ -423,8 +423,8 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - const value_type& operator()(const unsigned short i, - const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 6 + j]; } // end of operator() /*! @@ -432,7 +432,7 @@ namespace tfel::math { * In this case, the number of lines and columns * are deduced from the template parameter */ - RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct diff --git a/include/TFEL/Math/ST2toT2/T2toST2ST2toT2ProductExpr.hxx b/include/TFEL/Math/ST2toT2/T2toST2ST2toT2ProductExpr.hxx index 525f62bc7..b4a89626e 100644 --- a/include/TFEL/Math/ST2toT2/T2toST2ST2toT2ProductExpr.hxx +++ b/include/TFEL/Math/ST2toT2/T2toST2ST2toT2ProductExpr.hxx @@ -35,7 +35,7 @@ namespace tfel::math { template struct TFEL_VISIBILITY_LOCAL Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder< StensorDimeToSize()>::value * @@ -95,7 +95,7 @@ namespace tfel::math { template struct TFEL_VISIBILITY_LOCAL Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder< StensorDimeToSize()>::value * @@ -178,7 +178,7 @@ namespace tfel::math { template struct TFEL_VISIBILITY_LOCAL Expr> - : public ST2toST2Concept< + : public ST2toST2ConceptBase< Expr>>, public array_holder< StensorDimeToSize()>::value * diff --git a/include/TFEL/Math/ST2toT2/st2tot2.ixx b/include/TFEL/Math/ST2toT2/st2tot2.ixx index b7a8cffbf..245064d53 100644 --- a/include/TFEL/Math/ST2toT2/st2tot2.ixx +++ b/include/TFEL/Math/ST2toT2/st2tot2.ixx @@ -35,10 +35,9 @@ namespace tfel::math { } // end of st2tot2 template - template + template std::enable_if_t< - implementsST2toST2Concept() && - getSpaceDimension() == N && + getSpaceDimension() == N && getSpaceDimension() == N && isAssignableTo, numeric_type, @@ -52,33 +51,32 @@ namespace tfel::math { template template - std::enable_if_t() == N && - isAssignableTo, T>(), - Expr, StensorProductRightDerivativeExpr>> - st2tot2::tprd(const StensorType& a) { + TFEL_HOST_DEVICE constexpr std::enable_if_t< + getSpaceDimension() == N && + isAssignableTo, T>(), + Expr, StensorProductRightDerivativeExpr>> + st2tot2::tprd(const StensorType& a) noexcept { return Expr, StensorProductRightDerivativeExpr>(a); } template - template - std::enable_if_t< - implementsST2toST2Concept() && - getSpaceDimension() == N && + template + TFEL_HOST_DEVICE constexpr std::enable_if_t< + getSpaceDimension() == N && getSpaceDimension() == N && isAssignableTo, numeric_type, OpMult>, T>(), Expr, StensorProductRightDerivativeExpr>> - st2tot2::tprd(const StensorType& a, const ST2toST2Type& C) { + st2tot2::tprd(const StensorType& a, const ST2toST2Type& C) noexcept { return Expr, StensorProductRightDerivativeExpr>(a, C); } template - template - TFEL_MATH_INLINE2 void st2tot2::copy(const InputIterator src) { + TFEL_HOST_DEVICE constexpr void st2tot2::copy(const auto p) noexcept { tfel::fsalgo::copy::value * - StensorDimeToSize::value>::exe(src, *this); + StensorDimeToSize::value>::exe(p, *this); } } // end of namespace tfel::math diff --git a/include/TFEL/Math/Stensor/DecompositionInPositiveAndNegativeParts.hxx b/include/TFEL/Math/Stensor/DecompositionInPositiveAndNegativeParts.hxx index 8f60a63a3..683f3fc55 100644 --- a/include/TFEL/Math/Stensor/DecompositionInPositiveAndNegativeParts.hxx +++ b/include/TFEL/Math/Stensor/DecompositionInPositiveAndNegativeParts.hxx @@ -32,14 +32,15 @@ namespace tfel::math { * \pre PPType must implement the StensorConcept * \pre StensorType must implement the StensorConcept */ - template + template void computeStensorPositivePartAndDerivative( DPPType&, PPType&, const StensorType&, const numeric_type) // requires( - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == 1u && @@ -60,14 +61,15 @@ namespace tfel::math { * \pre PPType must implement the StensorConcept * \pre StensorType must implement the StensorConcept */ - template + template void computeStensorPositivePartAndDerivative( DPPType&, PPType&, const StensorType&, const numeric_type) // requires( - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == 2u && @@ -88,14 +90,15 @@ namespace tfel::math { * \pre PPType must implement the StensorConcept * \pre StensorType must implement the StensorConcept */ - template + template void computeStensorPositivePartAndDerivative( DPPType&, PPType&, const StensorType&, const numeric_type) // requires( - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == 3u && @@ -122,8 +125,8 @@ namespace tfel::math { * \pre NPType must implement the StensorConcept * \pre StensorType must implement the StensorConcept */ - template @@ -135,8 +138,6 @@ namespace tfel::math { const StensorType&, const numeric_type) // requires( - implementsST2toST2Concept() && - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && @@ -169,8 +170,8 @@ namespace tfel::math { * \pre NPType must implement the StensorConcept * \pre StensorType must implement the StensorConcept */ - template @@ -182,8 +183,6 @@ namespace tfel::math { const StensorType&, const numeric_type) // requires( - implementsST2toST2Concept() && - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && @@ -216,8 +215,8 @@ namespace tfel::math { * \pre NPType must implement the StensorConcept * \pre StensorType must implement the StensorConcept */ - template @@ -229,8 +228,6 @@ namespace tfel::math { const StensorType&, const numeric_type) // requires( - implementsST2toST2Concept() && - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && diff --git a/include/TFEL/Math/Stensor/DecompositionInPositiveAndNegativeParts.ixx b/include/TFEL/Math/Stensor/DecompositionInPositiveAndNegativeParts.ixx index ae1e0d457..69b3e9ea7 100644 --- a/include/TFEL/Math/Stensor/DecompositionInPositiveAndNegativeParts.ixx +++ b/include/TFEL/Math/Stensor/DecompositionInPositiveAndNegativeParts.ixx @@ -33,14 +33,15 @@ namespace tfel::math::internals { namespace tfel::math { - template + template void computeStensorPositivePartAndDerivative( DPPType& dpp, PPType& pp, const StensorType& s, const numeric_type eps) // requires( - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == 1u && @@ -85,14 +86,15 @@ namespace tfel::math { } } - template + template void computeStensorPositivePartAndDerivative( DPPType& dpp, PPType& pp, const StensorType& s, const numeric_type eps) // requires( - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == 2u && @@ -159,14 +161,15 @@ namespace tfel::math { } } - template + template void computeStensorPositivePartAndDerivative( DPPType& dpp, PPType& pp, const StensorType& s, const numeric_type eps) // requires( - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == 3u && @@ -309,8 +312,8 @@ namespace tfel::math { } } - template @@ -322,8 +325,6 @@ namespace tfel::math { const StensorType& s, const numeric_type eps) // requires( - implementsST2toST2Concept() && - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && @@ -390,8 +391,8 @@ namespace tfel::math { } } // end of computeStensorDecompositionInPositiveAndNegativeParts - template @@ -403,8 +404,6 @@ namespace tfel::math { const StensorType& s, const numeric_type eps) // requires( - implementsST2toST2Concept() && - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && @@ -494,8 +493,8 @@ namespace tfel::math { } } // end of computeStensorDecompositionInPositiveAndNegativeParts - template @@ -507,8 +506,6 @@ namespace tfel::math { const StensorType& s, const numeric_type eps) // requires( - implementsST2toST2Concept() && - implementsST2toST2Concept() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && getSpaceDimension() == getSpaceDimension() && diff --git a/include/TFEL/Math/Stensor/Internals/StensorChangeBasis.hxx b/include/TFEL/Math/Stensor/Internals/StensorChangeBasis.hxx index 7c88223db..fd64c18f6 100644 --- a/include/TFEL/Math/Stensor/Internals/StensorChangeBasis.hxx +++ b/include/TFEL/Math/Stensor/Internals/StensorChangeBasis.hxx @@ -28,13 +28,16 @@ namespace tfel::math::internals { template <> struct StensorChangeBasis<1u> { template - static void exe(T*, const tfel::math::rotation_matrix&) {} + TFEL_HOST_DEVICE static constexpr void exe( + T* const, const tfel::math::rotation_matrix>&) noexcept {} }; template <> struct StensorChangeBasis<2u> { template - static void exe(T* s, const tfel::math::rotation_matrix& m) { + TFEL_HOST_DEVICE static constexpr void exe( + T* const s, + const tfel::math::rotation_matrix>& m) noexcept { constexpr auto cste = Cste::sqrt2; T tmp[3]; // Works begin @@ -53,7 +56,9 @@ namespace tfel::math::internals { template <> struct StensorChangeBasis<3u> { template - static void exe(T* s, const tfel::math::rotation_matrix& m) { + TFEL_HOST_DEVICE static constexpr void exe( + T* const s, + const tfel::math::rotation_matrix>& m) noexcept { constexpr auto cste = Cste::sqrt2; T tmp[6]; const auto a = m(0, 0); diff --git a/include/TFEL/Math/Stensor/Internals/StensorComputeEigenVectorsDerivatives.hxx b/include/TFEL/Math/Stensor/Internals/StensorComputeEigenVectorsDerivatives.hxx index 495d0498e..6469356b3 100644 --- a/include/TFEL/Math/Stensor/Internals/StensorComputeEigenVectorsDerivatives.hxx +++ b/include/TFEL/Math/Stensor/Internals/StensorComputeEigenVectorsDerivatives.hxx @@ -41,21 +41,20 @@ namespace tfel::math::internals { template <> struct StensorComputeEigenTensorsDerivatives<1u> { - template - static TFEL_MATH_INLINE2 typename std::enable_if< - ((implementsST2toST2Concept()) && - (getSpaceDimension() == 1u) && - (isAssignableTo, - NumType, - OpDiv>::Result, - numeric_type>())), - void>::type - exe(ST2toST2Type& dn0_ds, + template + static TFEL_MATH_INLINE2 void exe( + ST2toST2Type& dn0_ds, ST2toST2Type& dn1_ds, ST2toST2Type& dn2_ds, const tfel::math::tvector<3u, NumType>&, const tfel::math::rotation_matrix&, - const NumType) { + const NumType) // + requires( + (getSpaceDimension() == 1u) && + (isAssignableTo, + NumType, + OpDiv>::Result, + numeric_type>())) { using namespace tfel::math; using namespace tfel::typetraits; typedef typename ComputeBinaryResult, NumType, @@ -69,21 +68,20 @@ namespace tfel::math::internals { template <> struct StensorComputeEigenTensorsDerivatives<2u> : public StensorComputeEigenTensorsDerivativesBase { - template - static TFEL_MATH_INLINE2 typename std::enable_if< - ((implementsST2toST2Concept()) && - (getSpaceDimension() == 2u) && - (isAssignableTo, - NumType, - OpDiv>::Result, - numeric_type>())), - void>::type - exe(ST2toST2Type& dn0_ds, + template + static TFEL_MATH_INLINE2 void exe( + ST2toST2Type& dn0_ds, ST2toST2Type& dn1_ds, ST2toST2Type& dn2_ds, const tfel::math::tvector<3u, NumType>& vp, const tfel::math::rotation_matrix& m, - const NumType eps) { + const NumType eps) // + requires( + (getSpaceDimension() == 2u) && + (isAssignableTo, + NumType, + OpDiv>::Result, + numeric_type>())) { typedef base_type base; typedef BinaryOperationResult InvNumType; constexpr auto icste = Cste::isqrt2; @@ -101,21 +99,20 @@ namespace tfel::math::internals { template <> struct StensorComputeEigenTensorsDerivatives<3u> : public StensorComputeEigenTensorsDerivativesBase { - template - static TFEL_MATH_INLINE2 typename std::enable_if< - ((implementsST2toST2Concept()) && - (getSpaceDimension() == 3u) && - (isAssignableTo, - NumType, - OpDiv>::Result, - numeric_type>())), - void>::type - exe(ST2toST2Type& dn0_ds, + template + static TFEL_MATH_INLINE2 void exe( + ST2toST2Type& dn0_ds, ST2toST2Type& dn1_ds, ST2toST2Type& dn2_ds, const tfel::math::tvector<3u, NumType>& vp, const tfel::math::rotation_matrix& m, - const NumType eps) { + const NumType eps) // + requires( + (getSpaceDimension() == 3u) && + (isAssignableTo, + NumType, + OpDiv>::Result, + numeric_type>())) { using namespace tfel::math; typedef base_type base; constexpr auto cste = Cste::isqrt2; diff --git a/include/TFEL/Math/Stensor/Internals/StensorComputeIsotropicFunctionDerivative.hxx b/include/TFEL/Math/Stensor/Internals/StensorComputeIsotropicFunctionDerivative.hxx index 94db96ae0..85bb90c58 100644 --- a/include/TFEL/Math/Stensor/Internals/StensorComputeIsotropicFunctionDerivative.hxx +++ b/include/TFEL/Math/Stensor/Internals/StensorComputeIsotropicFunctionDerivative.hxx @@ -31,20 +31,17 @@ namespace tfel::math::internals { * \param[in] eps: criterion value used to judge if two eigenvalues are * equals */ - template - static - typename std::enable_if<(implementsST2toST2Concept()) && - (getSpaceDimension() == 1u), - void>::type - exe2(ST2toST2Type& d, - const Function&, - const FunctionDerivative& df, - const tvector<3u, T>& vp, - const rotation_matrix&, - const T) { + static void exe2(ST2toST2Type& d, + const Function&, + const FunctionDerivative& df, + const tvector<3u, T>& vp, + const rotation_matrix&, + const T) // + requires(getSpaceDimension() == 1u) { using real = numeric_type; constexpr auto zero = real(0); d(0, 0) = df(vp[0]); @@ -64,17 +61,17 @@ namespace tfel::math::internals { * \param[in] eps: criterion value used to judge if two eigenvalues are * equals */ - template - static - typename std::enable_if<(implementsST2toST2Concept()) && - (getSpaceDimension() == 1u), - void>::type - exe(ST2toST2Type& d, - const tvector<3u, T1>&, - const tvector<3u, T2>& df, - const tvector<3u, T3>&, - const tmatrix<3u, 3u, base_type>&, - const T3) { + template + static void exe(ST2toST2Type& d, + const tvector<3u, T1>&, + const tvector<3u, T2>& df, + const tvector<3u, T3>&, + const tmatrix<3u, 3u, base_type>&, + const T3) // + requires(getSpaceDimension() == 1u) { using real = numeric_type; constexpr auto zero = real(0); d(0, 0) = df[0]; @@ -98,17 +95,17 @@ namespace tfel::math::internals { * \param[in] eps: criterion value used to judge if two eigenvalues are * equals */ - template - static - typename std::enable_if<(implementsST2toST2Concept()) && - (getSpaceDimension() == 2u), - void>::type - exe(ST2toST2Type& d, - const tvector<3u, T1>& f, - const tvector<3u, T2>& df, - const tvector<3u, T3>& vp, - const tmatrix<3u, 3u, base_type>& m, - const T3 eps) { + template + static void exe(ST2toST2Type& d, + const tvector<3u, T1>& f, + const tvector<3u, T2>& df, + const tvector<3u, T3>& vp, + const tmatrix<3u, 3u, base_type>& m, + const T3 eps) // + requires(getSpaceDimension() == 2u) { using real = numeric_type; using base = base_type; using tvector = tfel::math::tvector<3u, real>; @@ -138,20 +135,17 @@ namespace tfel::math::internals { * \param[in] eps: criterion value used to judge if two eigenvalues are * equals */ - template - static - typename std::enable_if<(implementsST2toST2Concept()) && - (getSpaceDimension() == 2u), - void>::type - exe2(ST2toST2Type& d, - const Function& f, - const FunctionDerivative& df, - const tvector<3u, T>& vp, - const rotation_matrix& m, - const T eps) { + static void exe2(ST2toST2Type& d, + const Function& f, + const FunctionDerivative& df, + const tvector<3u, T>& vp, + const rotation_matrix& m, + const T eps) // + requires(getSpaceDimension() == 2u) { const auto fv = map(f, vp); const auto dfv = map(df, vp); StensorComputeIsotropicFunctionDerivative::exe(d, fv, dfv, vp, m, eps); @@ -170,17 +164,17 @@ namespace tfel::math::internals { * \param[in] eps: criterion value used to judge if two eigenvalues are * equals */ - template - static - typename std::enable_if<(implementsST2toST2Concept()) && - (getSpaceDimension() == 3u), - void>::type - exe(ST2toST2Type& d, - const tvector<3u, T1>& f, - const tvector<3u, T2>& df, - const tvector<3u, T3>& vp, - const tmatrix<3u, 3u, base_type>& m, - const T3 eps) { + template + static void exe(ST2toST2Type& d, + const tvector<3u, T1>& f, + const tvector<3u, T2>& df, + const tvector<3u, T3>& vp, + const tmatrix<3u, 3u, base_type>& m, + const T3 eps) // + requires(getSpaceDimension() == 3u) { using real = numeric_type; using base = base_type; using tvector = tfel::math::tvector<3u, real>; @@ -258,20 +252,17 @@ namespace tfel::math::internals { * \param[in] eps: criterion value used to judge if two eigenvalues are * equals */ - template - static - typename std::enable_if<(implementsST2toST2Concept()) && - (getSpaceDimension() == 3u), - void>::type - exe2(ST2toST2Type& d, - const Function& f, - const FunctionDerivative& df, - const tvector<3u, T>& vp, - const rotation_matrix& m, - const T eps) { + static void exe2(ST2toST2Type& d, + const Function& f, + const FunctionDerivative& df, + const tvector<3u, T>& vp, + const rotation_matrix& m, + const T eps) // + requires(getSpaceDimension() == 3u) { const auto fv = map(f, vp); const auto dfv = map(df, vp); StensorComputeIsotropicFunctionDerivative::exe(d, fv, dfv, vp, m, eps); diff --git a/include/TFEL/Math/Stensor/StensorConcept.hxx b/include/TFEL/Math/Stensor/StensorConcept.hxx index b2aeda9e0..09acb325d 100644 --- a/include/TFEL/Math/Stensor/StensorConcept.hxx +++ b/include/TFEL/Math/Stensor/StensorConcept.hxx @@ -58,9 +58,7 @@ namespace tfel::math { (requires(const T t, const unsigned short i) { t[i]; }) && // (requires(const T t, const unsigned short i) { t(i); }); - /*! - * \brief partial specialisation for symmetric tensors - */ + //! \brief partial specialisation for symmetric tensors template struct ConceptRebind { //! \brief a simple alias diff --git a/include/TFEL/Math/Stensor/stensor.ixx b/include/TFEL/Math/Stensor/stensor.ixx index 5fa7045d7..8d4a1020e 100644 --- a/include/TFEL/Math/Stensor/stensor.ixx +++ b/include/TFEL/Math/Stensor/stensor.ixx @@ -262,10 +262,9 @@ namespace tfel::math { } // end of stensor::computeEigenTensors template - template + template TFEL_HOST_DEVICE std::enable_if_t< - (implementsST2toST2Concept()) && - (getSpaceDimension() == N) && + (getSpaceDimension() == N) && (isAssignableTo, T, OpDiv>, numeric_type>()), void> @@ -311,12 +310,11 @@ namespace tfel::math { } // end of stensor::computeIsotropicFunctionDerivative template - template TFEL_HOST_DEVICE std::enable_if_t< - (implementsST2toST2Concept()) && - (getSpaceDimension() == N) && + (getSpaceDimension() == N) && (isAssignableTo, T, OpDiv>, numeric_type>()), void> @@ -345,10 +343,9 @@ namespace tfel::math { } // end of stensor::computeIsotropicFunctionDerivative template - template + template TFEL_HOST_DEVICE std::enable_if_t< - (implementsST2toST2Concept()) && - (getSpaceDimension() == N) && + (getSpaceDimension() == N) && (isAssignableTo, T, OpDiv>, numeric_type>()), void> @@ -420,16 +417,14 @@ namespace tfel::math { return r; } // end of stensor::computeIsotropicFunctionAndDerivative - // ChangeBasis template - TFEL_HOST_DEVICE void stensor::changeBasis( - const rotation_matrix& m) { + TFEL_HOST_DEVICE constexpr void stensor::changeBasis( + const rotation_matrix>& m) noexcept { return tfel::math::internals::StensorChangeBasis::exe(this->v, m); } - // Return Id template - constexpr stensor> stensor::Id() { + constexpr stensor> stensor::Id() noexcept { static_assert((N == 1) || (N == 2) || (N == 3), "invalid space dimension"); constexpr auto zero = base_type{0}; constexpr auto one = base_type{1}; diff --git a/include/TFEL/Math/T2toST2/ST2toST2T2toST2ProductExpr.hxx b/include/TFEL/Math/T2toST2/ST2toST2T2toST2ProductExpr.hxx index eb9814bf1..e809b5fc0 100644 --- a/include/TFEL/Math/T2toST2/ST2toST2T2toST2ProductExpr.hxx +++ b/include/TFEL/Math/T2toST2/ST2toST2T2toST2ProductExpr.hxx @@ -22,15 +22,15 @@ namespace tfel::math { - //! Empty structure used for partial specialisation of the - //! Expr class + /*! + * \brief empty structure used for partial specialisation of the + * Expr class + */ template struct TFEL_VISIBILITY_LOCAL ST2toST2T2toST2ProductExpr { }; // end of struct ST2toST2T2toST2ProductExpr - /*! - * Partial specialisation - */ + //! \brief partial specialisation in 1D template struct TFEL_VISIBILITY_LOCAL Expr> @@ -49,9 +49,9 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const ST2toST2Type& a, const T2toST2Type2& b) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a, + const T2toST2Type2& b) noexcept { static_assert(implementsT2toST2Concept()); static_assert(getSpaceDimension() == 1u); static_assert(getSpaceDimension() == 1u); @@ -70,21 +70,17 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i, const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 3 + j]; } // end of operator() - /*! - * \return the runtime properties - */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + //! \return the runtime properties + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr - /*! - * Partial specialisation - */ + //! \brief partial specialisation in 2D template struct TFEL_VISIBILITY_LOCAL Expr> @@ -103,9 +99,9 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const ST2toST2Type& a, const T2toST2Type2& b) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a, + const T2toST2Type2& b) noexcept { static_assert(implementsT2toST2Concept()); static_assert(getSpaceDimension() == 2u); static_assert(getSpaceDimension() == 2u); @@ -155,21 +151,17 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i, const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 5 + j]; } // end of operator() - /*! - * \return the runtime properties - */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + //! \return the runtime properties + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr - /*! - * Partial specialisation - */ + //! \brief partial specialisation in 3D template struct TFEL_VISIBILITY_LOCAL Expr> @@ -188,9 +180,9 @@ namespace tfel::math { * \param[in] a : first term of the product * \param[in] b : second term of the product */ - template - TFEL_MATH_INLINE Expr(const ST2toST2Type& a, const T2toST2Type2& b) { - static_assert(implementsST2toST2Concept()); + template + TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a, + const T2toST2Type2& b) noexcept { static_assert(implementsT2toST2Concept()); static_assert(getSpaceDimension() == 3u); static_assert(getSpaceDimension() == 3u); @@ -308,14 +300,14 @@ namespace tfel::math { * \param[in] i : line index * \param[in] j : column index */ - TFEL_MATH_INLINE const value_type& operator()( - const unsigned short i, const unsigned short j) const { + TFEL_HOST_DEVICE constexpr const value_type& operator()( + const unsigned short i, const unsigned short j) const noexcept { return this->v[i * 9 + j]; } // end of operator() /*! * \return the runtime properties */ - TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const { + TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept { return RunTimeProperties(); } }; // end of struct Expr diff --git a/include/TFEL/Math/Tensor/tensor.ixx b/include/TFEL/Math/Tensor/tensor.ixx index d75397a3a..7c4b80a63 100644 --- a/include/TFEL/Math/Tensor/tensor.ixx +++ b/include/TFEL/Math/Tensor/tensor.ixx @@ -49,14 +49,14 @@ namespace tfel::math { } // end of void tensor::buildFromFortranMatrix template - constexpr ValueType tensor::operator()( - const typename tensor::size_type i) const { + TFEL_HOST_DEVICE constexpr ValueType tensor::operator()( + const typename tensor::size_type i) const noexcept { return GenericFixedSizeArrayBase::operator()(i); } // end of operator() template - constexpr ValueType& tensor::operator()( - const typename tensor::size_type i) { + TFEL_HOST_DEVICE constexpr ValueType& tensor::operator()( + const typename tensor::size_type i) noexcept { return GenericFixedSizeArrayBase::operator()(i); } @@ -98,26 +98,29 @@ namespace tfel::math { } // end of operator() template - TFEL_HOST_DEVICE void tensor::import(const base_type* const src) { + TFEL_HOST_DEVICE constexpr void tensor::import( + const base_type* const src) noexcept { tfel::fsalgo::transform::value>::exe( src, this->begin(), [](const auto& v) { return T(v); }); } template - void tensor::write(base_type* const t) const { + TFEL_HOST_DEVICE constexpr void tensor::write( + base_type* const t) const noexcept { tfel::fsalgo::transform::value>::exe( this->cbegin(), t, [](const auto& v) { return base_type_cast(v); }); } template - TFEL_HOST_DEVICE void tensor::changeBasis( + TFEL_HOST_DEVICE constexpr void tensor::changeBasis( const rotation_matrix& m) noexcept { const auto rt = change_basis(*this, m); *this = rt; } // end of changeBasis template - constexpr tensor> tensor::Id() { + TFEL_HOST_DEVICE constexpr tensor> + tensor::Id() noexcept { static_assert((N == 1) || (N == 2) || (N == 3), "invalid space dimension"); constexpr auto zero = base_type{0}; constexpr auto one = base_type{1}; @@ -130,10 +133,9 @@ namespace tfel::math { } // end of tensor::Id template - template - TFEL_HOST_DEVICE void tensor::copy(const InputIterator src) { - tfel::fsalgo::copy::value>::exe(src, *this); - } + TFEL_HOST_DEVICE constexpr void tensor::copy(const auto p) noexcept { + tfel::fsalgo::copy::value>::exe(p, *this); + } // end of copy template TFEL_HOST_DEVICE TFEL_MATH_INLINE2 std::enable_if_t(), void> diff --git a/include/TFEL/Math/st2tost2.hxx b/include/TFEL/Math/st2tost2.hxx index f35aba8d0..8932c8397 100644 --- a/include/TFEL/Math/st2tost2.hxx +++ b/include/TFEL/Math/st2tost2.hxx @@ -50,13 +50,11 @@ namespace tfel::math { * \brief partial specialisation of the `DerivativeTypeDispatcher` * metafunction. */ - template + template struct DerivativeTypeDispatcher { - static_assert(implementsST2toST2Concept(), - "template argument ST2toST2Type is not a st2tost2"); static_assert(isScalar(), "template argument ScalarType is not a scalar"); static_assert(isScalar>(), @@ -70,13 +68,11 @@ namespace tfel::math { * \brief partial specialisation of the `DerivativeTypeDispatcher` * metafunction. */ - template + template struct DerivativeTypeDispatcher { - static_assert(implementsST2toST2Concept(), - "template argument ST2toST2Type is not a st2tost2"); static_assert(isScalar(), "template argument ScalarType is not a scalar"); static_assert(isScalar>(), @@ -111,7 +107,7 @@ namespace tfel::math { template struct st2tost2 - : ST2toST2Concept>, + : ST2toST2ConceptBase>, GenericFixedSizeArray< st2tost2, FixedSizeRowMajorMatrixPolicy::value, @@ -129,27 +125,24 @@ namespace tfel::math { * \return the derivative of the square of a symmetric tensor */ template - static TFEL_MATH_INLINE std::enable_if_t< + TFEL_HOST_DEVICE static constexpr auto + dsquare(const StensorType&) noexcept requires( getSpaceDimension() == N && - isAssignableTo, ValueType>(), - Expr, StensorSquareDerivativeExpr>> - dsquare(const StensorType&); + isAssignableTo, ValueType>()); /*! * \param[in] s : tensor squared * \param[in] C : derivative of s * \return the derivative of the square of a symmetric tensor */ - template - static TFEL_MATH_INLINE std::enable_if_t< - implementsST2toST2Concept() && - getSpaceDimension() == N && - getSpaceDimension() == N && - isAssignableTo, - numeric_type, - OpMult>, - ValueType>(), - Expr, StensorSquareDerivativeExpr>> - dsquare(const StensorType&, const ST2toST2Type&); + template + TFEL_HOST_DEVICE static constexpr auto + dsquare(const StensorType&, const ST2toST2Type&) noexcept requires( + getSpaceDimension() == N && + getSpaceDimension() == N && + isAssignableTo, + numeric_type, + OpMult>, + ValueType>()); /*! * convert a T2toST2 to a ST2toST2 * \param[in] src : T2toST2 to be converted @@ -181,11 +174,11 @@ namespace tfel::math { tfel::math::st2tost2> stpd(const StensorType&); // - static constexpr st2tost2 Id(); - static constexpr st2tost2 IxI(); - static constexpr st2tost2 K(); - static constexpr st2tost2 M(); - static constexpr st2tost2 J(); + static constexpr st2tost2 Id() noexcept; + static constexpr st2tost2 IxI() noexcept; + static constexpr st2tost2 K() noexcept; + static constexpr st2tost2 M() noexcept; + static constexpr st2tost2 J() noexcept; // TFEL_MATH_FIXED_SIZE_ARRAY_DEFAULT_METHODS(st2tost2, GenericFixedSizeArrayBase); @@ -193,10 +186,13 @@ namespace tfel::math { using GenericFixedSizeArrayBase::operator[]; using GenericFixedSizeArrayBase::operator(); //! \brief import values from an external memory location - void import(const base_type* const); - - template - TFEL_MATH_INLINE2 void copy(const InputIterator src); + TFEL_HOST_DEVICE constexpr void import( + const base_type* const) noexcept; + /*! + * \brief copy the values from a range + * \param[in] p: iterator + */ + TFEL_HOST_DEVICE constexpr void copy(const auto) noexcept; }; /*! @@ -219,75 +215,46 @@ namespace tfel::math { * \param[in] s : st2tost2 * \param[in] r : rotation matrix */ - template - TFEL_MATH_INLINE2 std::enable_if_t< - implementsST2toST2Concept(), - st2tost2(), numeric_type>> + template + TFEL_HOST_DEVICE constexpr st2tost2(), + numeric_type> change_basis(const ST2toST2Type&, - const rotation_matrix>&); + const rotation_matrix>&) noexcept; /*! * \return the invert of a st2tost2 * \param[in] s : st2tost2 to be inverted */ - template - TFEL_MATH_INLINE2 std::enable_if_t< - implementsST2toST2Concept(), - st2tost2(), - BinaryOperationResult>, - numeric_type, - OpDiv>>> - invert(const ST2toST2Type&); + TFEL_HOST constexpr auto invert(const ST2toST2Concept auto&); /*! - * \return the push_forward of a st2tost2: + * \return the push-forward of a fourth order tensor: * \f[ * Ct_{ijkl}=F_{im}F_{jn}F_{kp}F_{lq}C_{mnpq} * \f] * \param[in] C: input * \param[in] F: deformation gradient */ - template - std::enable_if_t() && - getSpaceDimension() == - getSpaceDimension(), - st2tost2(), - BinaryOperationResult, - numeric_type, - OpMult>>> - push_forward(const ST2toST2Type&, const TensorType&); - - template - std::enable_if_t() && - getSpaceDimension() == - getSpaceDimension(), - st2tost2(), - BinaryOperationResult, - numeric_type, - OpMult>>> - pull_back(const ST2toST2Type&, const TensorType&); + template + TFEL_HOST_DEVICE constexpr auto push_forward( + const ST2toST2Type&, + const TensorType&) noexcept requires(getSpaceDimension() == + getSpaceDimension()); /*! - * \brief compute the second derivative of determinant of the - * deviator of a symmetric tensor with respect to this tensor. - * - * Let \f$\underline{s}\f$ be a symmetric tensor and \f$J_{3}\f$ - * be the determinant of \f$\underline{s}'\f$ the deviator of - * \f$\underline{s}\f$: - * \f[ - * J_{3} = \mathrm{det}\left(\underline{s}'\right) - * = - * \mathrm{det}\left(\underline{s}-\mathrm{tr}\left(\underline{s}'\right)\,\underline{I}\right) - * \f] - * - * This function computes \f$\displaystyle\frac{\partial^{2} J_{3}}{\partial - * \underline{\sigma}^{2}}\f$. - * - * \[ \param[in] s: tensor + * \return the pull-back of a fourth order tensor + * \param[in] C: input + * \param[in] F: deformation gradient + */ + template + TFEL_HOST constexpr auto pull_back( + const ST2toST2Type&, + const TensorType&) requires(getSpaceDimension() == + getSpaceDimension()); + /*! + * \brief compute the second derivative of determinant of a + * \param[in] s: tensor */ - template - std::enable_if_t< - isScalar>(), - st2tost2(), numeric_type>> - computeDeviatorDeterminantSecondDerivative(const StensorType&); + TFEL_HOST_DEVICE constexpr auto computeDeterminantSecondDerivative( + const StensorConcept auto&) noexcept; /*! * \brief compute the second derivative of the determinant of the * deviator of symmetric tensor. @@ -302,15 +269,12 @@ namespace tfel::math { * \f] * * This function computes \f$\displaystyle\frac{\partial^{2} J_{3}}{\partial - \underline{\sigma}^{2}}\f$. + * \underline{\sigma}^{2}}\f$. * * \param[in] s: tensor */ - template - std::enable_if_t< - isScalar>(), - st2tost2(), numeric_type>> - computeDeviatorDeterminantSecondDerivative(const StensorType&); + TFEL_HOST_DEVICE constexpr auto computeDeviatorDeterminantSecondDerivative( + const StensorConcept auto&) noexcept; } // end of namespace tfel::math diff --git a/include/TFEL/Math/st2tot2.hxx b/include/TFEL/Math/st2tot2.hxx index d9381f04e..d0deb0937 100644 --- a/include/TFEL/Math/st2tot2.hxx +++ b/include/TFEL/Math/st2tot2.hxx @@ -128,10 +128,9 @@ namespace tfel::math { * \param[in] C : derivative of the first tensor * \return the left part of the derivative of a tensor product */ - template + template static TFEL_MATH_INLINE std::enable_if_t< - implementsST2toST2Concept() && - getSpaceDimension() == N && + getSpaceDimension() == N && getSpaceDimension() == N && isAssignableTo, numeric_type, @@ -144,27 +143,26 @@ namespace tfel::math { * \return the right part of the derivative of a tensor product */ template - static TFEL_MATH_INLINE std::enable_if_t< + TFEL_HOST_DEVICE static constexpr std::enable_if_t< getSpaceDimension() == N && isAssignableTo, ValueType>(), Expr, StensorProductRightDerivativeExpr>> - tprd(const StensorType&); + tprd(const StensorType&) noexcept; /*! * \param[in] A : first tensor of the product * \param[in] C : derivative of the first tensor * \return the right part of the derivative of a tensor product */ - template - static TFEL_MATH_INLINE std::enable_if_t< - implementsST2toST2Concept() && - getSpaceDimension() == N && + template + TFEL_HOST_DEVICE static constexpr std::enable_if_t< + getSpaceDimension() == N && getSpaceDimension() == N && isAssignableTo, numeric_type, OpMult>, ValueType>(), Expr, StensorProductRightDerivativeExpr>> - tprd(const StensorType&, const ST2toST2Type&); + tprd(const StensorType&, const ST2toST2Type&) noexcept; // TFEL_MATH_FIXED_SIZE_ARRAY_DEFAULT_METHODS(st2tot2, GenericFixedSizeArrayBase); @@ -174,8 +172,7 @@ namespace tfel::math { //! \brief import values from an external memory location void import(const base_type* const); // - template - TFEL_MATH_INLINE2 void copy(const InputIterator src); + TFEL_HOST_DEVICE constexpr void copy(const auto) noexcept; }; /*! diff --git a/include/TFEL/Math/stensor.hxx b/include/TFEL/Math/stensor.hxx index 235da43c7..01010f0e5 100644 --- a/include/TFEL/Math/stensor.hxx +++ b/include/TFEL/Math/stensor.hxx @@ -278,12 +278,13 @@ namespace tfel::math { VectorType&, const ValueType) const; /*! * \brief change basis + * \param[in] m: rotation matrix */ - TFEL_HOST_DEVICE TFEL_MATH_INLINE2 void changeBasis( - const rotation_matrix&); + TFEL_HOST_DEVICE constexpr void changeBasis( + const rotation_matrix>&) noexcept; //! \return the identity - TFEL_HOST_DEVICE - TFEL_MATH_INLINE static constexpr stensor> Id(); + TFEL_HOST_DEVICE static constexpr stensor> + Id() noexcept; /*! * copy the value from a container */ @@ -517,10 +518,9 @@ namespace tfel::math { * \param[in] m: eigen vectors * \param[in] eps: numerical parameter for regularisation */ - template + template TFEL_HOST_DEVICE static std::enable_if_t< - (implementsST2toST2Concept()) && - (getSpaceDimension() == N) && + (getSpaceDimension() == N) && (isAssignableTo< BinaryOperationResult, ValueType, OpDiv>, numeric_type>()), @@ -577,10 +577,9 @@ namespace tfel::math { * \param[in] eps: criterion value used to judge if two eigenvalues are * equals */ - template + template TFEL_HOST_DEVICE static std::enable_if_t< - (implementsST2toST2Concept()) && - (getSpaceDimension() == N) && + (getSpaceDimension() == N) && (isAssignableTo< BinaryOperationResult, ValueType, OpDiv>, numeric_type>()), @@ -620,12 +619,11 @@ namespace tfel::math { * \param[in] eps: criterion value used to judge if two eigenvalues are * equals */ - template TFEL_HOST_DEVICE static std::enable_if_t< - (implementsST2toST2Concept()) && - (getSpaceDimension() == N) && + (getSpaceDimension() == N) && (isAssignableTo< BinaryOperationResult, ValueType, OpDiv>, numeric_type>()), diff --git a/include/TFEL/Math/tensor.hxx b/include/TFEL/Math/tensor.hxx index f71b19240..25f0fc30a 100644 --- a/include/TFEL/Math/tensor.hxx +++ b/include/TFEL/Math/tensor.hxx @@ -84,7 +84,8 @@ namespace tfel::math { /*! * \return the identity tensor */ - TFEL_HOST_DEVICE static constexpr tensor> Id(); + TFEL_HOST_DEVICE static constexpr tensor> + Id() noexcept; /*! * \brief Build a tensor from a fortran matrix. * \param[in] t: tensor to be filled @@ -114,13 +115,13 @@ namespace tfel::math { * \param[in] i: index */ TFEL_HOST_DEVICE constexpr ValueType operator()( - const typename tensor::size_type) const; + const typename tensor::size_type) const noexcept; /*! * \brief access operator * \param[in] i: index */ TFEL_HOST_DEVICE constexpr ValueType& operator()( - const typename tensor::size_type); + const typename tensor::size_type) noexcept; /*! * \brief matrix-like access operator * \param[in] i: row number @@ -130,15 +131,16 @@ namespace tfel::math { const typename tensor::size_type, const typename tensor::size_type) const; //! \brief write to an external memory location - TFEL_MATH_INLINE2 void write(base_type* const) const; + TFEL_HOST_DEVICE constexpr void write(base_type* const) const + noexcept; //! \brief import values from an external memory location - TFEL_HOST_DEVICE void import(const base_type* const); + TFEL_HOST_DEVICE constexpr void import( + const base_type* const) noexcept; //! \brief change basis - TFEL_HOST_DEVICE TFEL_MATH_INLINE2 void changeBasis( + TFEL_HOST_DEVICE constexpr void changeBasis( const rotation_matrix&) noexcept; - - template - TFEL_HOST_DEVICE TFEL_MATH_INLINE2 void copy(const InputIterator src); + //! \brief copy from a range + TFEL_HOST_DEVICE constexpr void copy(const auto) noexcept; }; // end of struct tensor /*!