diff --git a/src/common/T4_map_proxy.hh b/src/common/T4_map_proxy.hh index c01c71c..596130d 100644 --- a/src/common/T4_map_proxy.hh +++ b/src/common/T4_map_proxy.hh @@ -1,112 +1,152 @@ /** * file T4_map_proxy.hh * * @author Till Junge * * @date 19 Nov 2017 * * @brief Map type to allow fourth-order tensor-like maps on 2D matrices * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include +#include "Eigen/src/Core/util/Constants.h" #ifndef T4_MAP_PROXY_H #define T4_MAP_PROXY_H namespace muSpectre { /* ---------------------------------------------------------------------- */ /** Proxy class mapping a fourth-order tensor onto a 2D matrix (in - order to avoid the use of Eigen::Tensor + order to avoid the use of Eigen::Tensor. This class is, however + byte-compatible with Tensors (i.e., you can map this onto a + tensor instead of a matrix) **/ template + int MapOptions=Eigen::Unaligned, + typename StrideType=Eigen::Stride<0, 0>> class T4Map: - public Eigen::MapBase, MapOptions> { - { + public Eigen::MapBase> { public: - typedef MapBase Base; - EIGEN_DENSE_PUBLIC_INTERFACE(T4Map); - - typedef typename Base::PointerType PointerType; - typedef PointerType PointerArgType; - EIGEN_DEVICE_FUNC - inline PointerType cast_to_pointer_type(PointerArgType ptr) { return ptr; } - - EIGEN_DEVICE_FUNC - inline Index innerStride() const - { - return StrideType::InnerStrideAtCompileTime != 0 ? m_stride.inner() : 1; - } - - EIGEN_DEVICE_FUNC - inline Index outerStride() const - { - return StrideType::OuterStrideAtCompileTime != 0 ? m_stride.outer() - : IsVectorAtCompileTime ? this->size() - : int(Flags)&RowMajorBit ? this->cols() - : this->rows(); - } - - /** Constructor in the fixed-size case. - * - * \param dataPtr pointer to the array to map - * \param stride optional Stride object, passing the strides. - */ - EIGEN_DEVICE_FUNC - explicit inline T4Map(PointerArgType dataPtr, const StrideType& stride = StrideType()) - : Base(cast_to_pointer_type(dataPtr)), m_stride(stride) - { - PlainObjectType::Base::_check_template_params(); - } - - EIGEN_INHERIT_ASSIGNMENT_OPERATORS(T4Map); - - /** My accessor to mimick tensorial access - **/ - inline const Scalar& operator()(Dim_t i, Dim_t j, Dim_t k, Dim_t l ) const { - constexpr auto myColStride{ - (colStride() == 1) ? colStride(): colStride()/Dim}; - constexpr auto myRowStride{ - (rowStride() == 1) ? rowStride(): rowStride()/Dim}; - return this->operator()(i * myRowStride + j * myColStride, - k * myRowStride + l * myColStride); - } - - inline Scalar& operator()(Dim_t i, Dim_t j, Dim_t k, Dim_t l ) { - constexpr auto myColStride{ - (colStride() == 1) ? colStride(): colStride()/Dim}; - constexpr auto myRowStride{ - (rowStride() == 1) ? rowStride(): rowStride()/Dim}; - return this->operator()(i * myRowStride + j * myColStride, - k * myRowStride + l * myColStride); - } - - protected: - StrideType m_stride; - }; + typedef Eigen::MapBase Base; + EIGEN_DENSE_PUBLIC_INTERFACE(T4Map); + + using PlainObjectType = Eigen::Matrix; + using Base::colStride; + using Base::rowStride; + using Base::IsRowMajor; + typedef typename Base::PointerType PointerType; + typedef PointerType PointerArgType; + EIGEN_DEVICE_FUNC + inline PointerType cast_to_pointer_type(PointerArgType ptr) { return ptr; } + + EIGEN_DEVICE_FUNC + inline Eigen::Index innerStride() const + { + return StrideType::InnerStrideAtCompileTime != 0 ? m_stride.inner() : 1; + } + + EIGEN_DEVICE_FUNC + inline Eigen::Index outerStride() const + { + return StrideType::OuterStrideAtCompileTime != 0 ? m_stride.outer() + : IsVectorAtCompileTime ? this->size() + : int(Flags)&Eigen::RowMajorBit ? this->cols() + : this->rows(); + } + + /** Constructor in the fixed-size case. + * + * \param dataPtr pointer to the array to map + * \param stride optional Stride object, passing the strides. + */ + EIGEN_DEVICE_FUNC + explicit inline T4Map(PointerArgType dataPtr, const StrideType& stride = StrideType()) + : Base(cast_to_pointer_type(dataPtr)), m_stride(stride), + map(cast_to_pointer_type(dataPtr)) + { + PlainObjectType::Base::_check_template_params(); + } + + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(T4Map); + + /** My accessor to mimick tensorial access + **/ + inline const Scalar& operator()(Dim_t i, Dim_t j, Dim_t k, Dim_t l ) const { + constexpr auto myColStride{ + (colStride() == 1) ? colStride(): colStride()/Dim}; + constexpr auto myRowStride{ + (rowStride() == 1) ? rowStride(): rowStride()/Dim}; + return this->operator()(i * myRowStride + j * myColStride, + k * myRowStride + l * myColStride); + } + + + inline Scalar& operator()(Dim_t i, Dim_t j, Dim_t k, Dim_t l ) { + const auto myColStride{ + (colStride() == 1) ? colStride(): colStride()/Dim}; + const auto myRowStride{ + (rowStride() == 1) ? rowStride(): rowStride()/Dim}; + return this->map(i * myRowStride + j * myColStride, + k * myRowStride + l * myColStride); + } + + protected: + StrideType m_stride; + Eigen::Map map; + }; } // muSpectre +namespace Eigen { + //! forward declarations + template struct traits; + + /* ---------------------------------------------------------------------- */ + namespace internal { + template + struct traits > + : public traits> + { + using PlainObjectType = Matrix; + typedef traits TraitsBase; + enum { + InnerStrideAtCompileTime = StrideType::InnerStrideAtCompileTime == 0 + ? int(PlainObjectType::InnerStrideAtCompileTime) + : int(StrideType::InnerStrideAtCompileTime), + OuterStrideAtCompileTime = StrideType::OuterStrideAtCompileTime == 0 + ? int(PlainObjectType::OuterStrideAtCompileTime) + : int(StrideType::OuterStrideAtCompileTime), + Alignment = int(MapOptions)&int(AlignedMask), + Flags0 = TraitsBase::Flags & (~NestByRefBit), + Flags = is_lvalue::value ? int(Flags0) : (int(Flags0) & ~LvalueBit) + }; + private: + enum { Options }; // Expressions don't have Options + }; + } // namespace internal +} // namespace Eigen + #endif /* T4_MAP_PROXY_H */ diff --git a/src/common/common.cc b/src/common/common.cc index b03485d..0746c45 100644 --- a/src/common/common.cc +++ b/src/common/common.cc @@ -1,142 +1,86 @@ /** * file common.cc * * @author Till Junge * * @date 15 Nov 2017 * * @brief Implementation for common functions * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include "common/common.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ std::ostream & operator<<(std::ostream & os, Formulation f) { switch (f) { case Formulation::small_strain: {os << "small_strain"; break;} case Formulation::finite_strain: {os << "finite_strain"; break;} default: throw std::runtime_error ("unknown formulation."); break; } return os; } /* ---------------------------------------------------------------------- */ std::ostream & operator<<(std::ostream & os, StressMeasure s) { switch (s) { case StressMeasure::Cauchy: {os << "Cauchy"; break;} case StressMeasure::PK1: {os << "PK1"; break;} case StressMeasure::PK2: {os << "PK2"; break;} case StressMeasure::Kirchhoff: {os << "Kirchhoff"; break;} case StressMeasure::Biot: {os << "Biot"; break;} case StressMeasure::Mandel: {os << "Mandel"; break;} default: throw std::runtime_error ("a stress measure must be missing"); break; } return os; } /* ---------------------------------------------------------------------- */ std::ostream & operator<<(std::ostream & os, StrainMeasure s) { switch (s) { case StrainMeasure::Gradient: {os << "Gradient"; break;} case StrainMeasure::Infinitesimal: {os << "Infinitesimal"; break;} case StrainMeasure::GreenLagrange: {os << "Green-Lagrange"; break;} case StrainMeasure::Biot: {os << "Biot"; break;} case StrainMeasure::Log: {os << "Logarithmic"; break;} case StrainMeasure::Almansi: {os << "Almansi"; break;} case StrainMeasure::RCauchyGreen: {os << "Right Cauchy-Green"; break;} case StrainMeasure::LCauchyGreen: {os << "Left Cauchy-Green"; break;} default: throw std::runtime_error ("a strain measure must be missing"); } return os; } - /* ---------------------------------------------------------------------- */ - /** Compile-time functions to set the stress and strain measures - stored by mu_spectre depending on the formulation - **/ - constexpr StrainMeasure get_stored_strain_type(Formulation form) { - switch (form) { - case Formulation::finite_strain: { - return StrainMeasure::Gradient; - break; - } - case Formulation::small_strain: { - return StrainMeasure::Infinitesimal; - break; - } - default: - return StrainMeasure::__nostrain__; - break; - } - } - - constexpr StressMeasure get_stored_stress_type(Formulation form) { - switch (form) { - case Formulation::finite_strain: { - return StressMeasure::PK1; - break; - } - case Formulation::small_strain: { - return StressMeasure::Cauchy; - break; - } - default: - return StressMeasure::__nostress__; - break; - } - } - - /* ---------------------------------------------------------------------- */ - constexpr StrainMeasure get_formulation_strain_type(Formulation form, - StrainMeasure expected) { - switch (form) { - case Formulation::finite_strain: { - return expected; - break; - } - case Formulation::small_strain: { - return get_stored_strain_type(form); - break; - } - default: - return StrainMeasure::__nostrain__; - break; - } - } - - - } // muSpectre diff --git a/src/common/common.hh b/src/common/common.hh index cac9d3f..90fac55 100644 --- a/src/common/common.hh +++ b/src/common/common.hh @@ -1,122 +1,166 @@ /** * file common.hh * * @author Till Junge * * @date 01 May 2017 * * @brief Small definitions of commonly used types througout µSpectre * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include #ifndef COMMON_H #define COMMON_H namespace muSpectre { //! Eigen uses signed integers for dimensions. For consistency, µSpectre uses //! them througout the code using Dim_t = int;// needs to represent -1 for eigen const Dim_t oneD{1}; const Dim_t twoD{2}; const Dim_t threeD{3}; + const Dim_t secondOrder{2}; + const Dim_t fourthOrder{4}; //! Ccoord_t are cell coordinates, i.e. integer coordinates template using Ccoord_t = std::array; template std::ostream & operator << (std::ostream & os, const Ccoord_t & index) { os << "("; for (size_t i = 0; i < dim-1; ++i) { os << index[i] << ", "; } os << index.back() << ")"; return os; } //! Scalar types used for mathematical calculations using Uint = unsigned int; using Int = int; using Real = double; using Complex = std::complex; //! compile-time potentiation required for field-size computations template constexpr I ipow(I base, I exponent) { static_assert(std::is_integral::value, "Type must be integer"); I retval{1}; for (I i = 0; i < exponent; ++i) { retval *= base; } return retval; } //! continuum mechanics flags enum class Formulation{finite_strain, small_strain}; std::ostream & operator<<(std::ostream & os, Formulation f); /* ---------------------------------------------------------------------- */ //! Material laws can declare which type of stress measure they provide, //! and µSpectre will handle conversions enum class StressMeasure { Cauchy, PK1, PK2, Kirchhoff, Biot, Mandel, __nostress__}; std::ostream & operator<<(std::ostream & os, StressMeasure s); /* ---------------------------------------------------------------------- */ //! Material laws can declare which type of strain measure they require and //! µSpectre will provide it enum class StrainMeasure { Gradient, Infinitesimal, GreenLagrange, Biot, Log, Almansi, RCauchyGreen, LCauchyGreen, __nostrain__}; std::ostream & operator<<(std::ostream & os, StrainMeasure s); /* ---------------------------------------------------------------------- */ /** Compile-time functions to set the stress and strain measures stored by mu_spectre depending on the formulation **/ - constexpr StrainMeasure get_stored_strain_type(Formulation form); - constexpr StressMeasure get_stored_stress_type(Formulation form); + constexpr StrainMeasure get_stored_strain_type(Formulation form) { + switch (form) { + case Formulation::finite_strain: { + return StrainMeasure::Gradient; + break; + } + case Formulation::small_strain: { + return StrainMeasure::Infinitesimal; + break; + } + default: + return StrainMeasure::__nostrain__; + break; + } + } + constexpr StressMeasure get_stored_stress_type(Formulation form) { + switch (form) { + case Formulation::finite_strain: { + return StressMeasure::PK1; + break; + } + case Formulation::small_strain: { + return StressMeasure::Cauchy; + break; + } + default: + return StressMeasure::__nostress__; + break; + } + } /* ---------------------------------------------------------------------- */ /** Compile-time functions to get the stress and strain measures after they may have been modified by choosing a formulation. For instance, a law that expecs a Green-Lagrange strain as input will get the infinitesimal strain tensor instead in a small strain computation **/ constexpr StrainMeasure get_formulation_strain_type(Formulation form, - StrainMeasure expected); + StrainMeasure expected) { + switch (form) { + case Formulation::finite_strain: { + return expected; + break; + } + case Formulation::small_strain: { + return get_stored_strain_type(form); + break; + } + default: + return StrainMeasure::__nostrain__; + break; + } + } } // muSpectre #ifndef EXPLICITLY_TURNED_ON_CXX17 #include "common/utilities.hh" #endif #endif /* COMMON_H */ diff --git a/src/common/field_map_matrixlike.hh b/src/common/field_map_matrixlike.hh index 2bd0a71..cad2979 100644 --- a/src/common/field_map_matrixlike.hh +++ b/src/common/field_map_matrixlike.hh @@ -1,195 +1,204 @@ /** * file field_map_matrixlike.hh * * @author Till Junge * * @date 26 Sep 2017 * * @brief Eigen-Matrix and -Array maps over strongly typed fields * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_MAP_MATRIXLIKE_H #define FIELD_MAP_MATRIXLIKE_H #include #include #include "common/field_map_base.hh" +#include "common/T4_map_proxy.hh" namespace muSpectre { namespace internal { /* ---------------------------------------------------------------------- */ enum class Map_t{Matrix, Array, T4Matrix}; template struct NameWrapper { const static std::string field_info_root; }; template<> const std::string NameWrapper::field_info_root{"Array"}; template<> const std::string NameWrapper::field_info_root{"Matrix"}; template<> const std::string NameWrapper::field_info_root{"T4Matrix"}; /* ---------------------------------------------------------------------- */ template class MatrixLikeFieldMap: public FieldMap { public: using Parent = FieldMap; using Ccoord = Ccoord_t; using value_type = EigenArray; using reference = value_type; // since it's a resource handle using const_reference = const value_type; using size_type = typename Parent::size_type; using pointer = std::unique_ptr; using TypedField = typename Parent::TypedField; using Field = typename TypedField::Parent; using iterator = typename Parent::template iterator; using const_iterator= typename Parent::template iterator; using reverse_iterator = std::reverse_iterator; using const_reverse_iterator = std::reverse_iterator; friend iterator; //! Default constructor MatrixLikeFieldMap() = delete; MatrixLikeFieldMap(Field & field); //! Copy constructor MatrixLikeFieldMap(const MatrixLikeFieldMap &other) = delete; //! Move constructor MatrixLikeFieldMap(MatrixLikeFieldMap &&other) noexcept = delete; //! Destructor virtual ~MatrixLikeFieldMap() noexcept = default; //! Copy assignment operator MatrixLikeFieldMap& operator=(const MatrixLikeFieldMap &other) = delete; //! Move assignment operator MatrixLikeFieldMap& operator=(MatrixLikeFieldMap &&other) noexcept = delete; //! give human-readable field map type inline std::string info_string() const override final; //! member access template inline ref operator[](size_type index); inline reference operator[](const Ccoord& ccoord); //! return an iterator to head of field for ranges inline iterator begin(){return iterator(*this);} inline const_iterator cbegin(){return const_iterator(*this);} //! return an iterator to tail of field for ranges inline iterator end(){return iterator(*this, false);}; inline const_iterator cend(){return const_iterator(*this, false);} protected: //! for sad, legacy iterator use inline pointer ptr_to_val_t(size_type index); const static std::string field_info_root; private: }; template const std::string MatrixLikeFieldMap:: field_info_root{NameWrapper::field_info_root}; /* ---------------------------------------------------------------------- */ template MatrixLikeFieldMap:: MatrixLikeFieldMap(Field & field) :Parent(field) { this->check_compatibility(); } /* ---------------------------------------------------------------------- */ //! human-readable field map type template std::string MatrixLikeFieldMap:: info_string() const { std::stringstream info; info << field_info_root << "(" << typeid(typename EigenArray::value_type).name() << ", " << EigenArray::RowsAtCompileTime << "x" << EigenArray::ColsAtCompileTime << ")"; return info.str(); } /* ---------------------------------------------------------------------- */ //! member access template template ref MatrixLikeFieldMap:: operator[](size_type index) { return ref(this->get_ptr_to_entry(index)); } template typename MatrixLikeFieldMap::reference MatrixLikeFieldMap:: operator[](const Ccoord & ccoord) { auto && index = this->collection.get_index(ccoord); return reference(this->get_ptr_to_entry(std::move(index))); } /* ---------------------------------------------------------------------- */ template typename MatrixLikeFieldMap::pointer MatrixLikeFieldMap:: ptr_to_val_t(size_type index) { return std::make_unique (this->get_ptr_to_entry(std::move(index))); } } // internal /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen matrix map as iterate template using MatrixFieldMap = internal::MatrixLikeFieldMap >, internal::Map_t::Matrix>; + /* ---------------------------------------------------------------------- */ + //! short-hand for an Eigen matrix map as iterate + template + using T4MatrixFieldMap = internal::MatrixLikeFieldMap + , + internal::Map_t::T4Matrix>; + /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen array map as iterate template using ArrayFieldMap = internal::MatrixLikeFieldMap >, internal::Map_t::Array>; } // muSpectre #endif /* FIELD_MAP_MATRIXLIKE_H */ diff --git a/src/common/tensor_algebra.hh b/src/common/tensor_algebra.hh index d798e3e..df3eda6 100644 --- a/src/common/tensor_algebra.hh +++ b/src/common/tensor_algebra.hh @@ -1,118 +1,125 @@ /** * file tensor_algebra.hh * * @author Till Junge * * @date 05 Nov 2017 * * @brief collection of compile-time quantities and algrebraic functions for * tensor operations * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef TENSOR_ALGEBRA_H #define TENSOR_ALGEBRA_H #include #include #include #include "common/common.hh" namespace muSpectre { namespace Tensors { template using Tens2_t = Eigen::TensorFixedSize>; template using Tens4_t = Eigen::TensorFixedSize>; //----------------------------------------------------------------------------// //! compile-time second-order identity template constexpr inline Tens2_t I2() { Tens2_t T; using Mat_t = Eigen::Matrix; Eigen::Map(&T(0, 0)) = Mat_t::Identity(); return T; } /* ---------------------------------------------------------------------- */ //! Check whether a given expression represents a Tensor specified order template struct is_tensor { constexpr static bool value = (std::is_convertible >::value || std::is_convertible >::value || std::is_convertible >::value); }; /* ---------------------------------------------------------------------- */ /** compile-time outer tensor product as defined by Curnier * R_ijkl = A_ij.B_klxx * 0123 01 23 */ template constexpr inline decltype(auto) outer(T1 && A, T2 && B) { // Just make sure that the right type of parameters have been given constexpr Dim_t order{2}; static_assert(is_tensor::value, "T1 needs to be convertible to a second order Tensor"); static_assert(is_tensor::value, "T2 needs to be convertible to a second order Tensor"); // actual function std::array, 0> dims{}; return A.contract(B, dims); } /* ---------------------------------------------------------------------- */ /** compile-time underlined outer tensor product as defined by Curnier * R_ijkl = A_ik.B_jlxx * 0123 02 13 * 0213 01 23 <- this defines the shuffle order */ template constexpr inline decltype(auto) outer_under(T1 && A, T2 && B) { constexpr size_t order{4}; return outer(A, B).shuffle(std::array{0, 2, 1, 3}); } /* ---------------------------------------------------------------------- */ /** compile-time overlined outer tensor product as defined by Curnier * R_ijkl = A_il.B_jkxx * 0123 03 12 * 0231 01 23 <- this defines the shuffle order */ template constexpr inline decltype(auto) outer_over(T1 && A, T2 && B) { constexpr size_t order{4}; return outer(A, B).shuffle(std::array{0, 2, 3, 1}); } + //! compile-time fourth-order symmetrising identity + template + constexpr inline Tens4_t I4S() { + auto I = I2(); + return 0.5*(outer_under(I, I) + outer_over(I, I)); + } + } // Tensors } // muSpectre #endif /* TENSOR_ALGEBRA_H */ diff --git a/src/common/utilities.hh b/src/common/utilities.hh index 70adc13..37b2f15 100644 --- a/src/common/utilities.hh +++ b/src/common/utilities.hh @@ -1,143 +1,143 @@ /** * file utilities.hh * * @author Till Junge * * @date 17 Nov 2017 * * @brief additions to the standard name space to anticipate C++17 features * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef UTILITIES_H #define UTILITIES_H namespace std { namespace detail { template struct is_reference_wrapper : std::false_type {}; template struct is_reference_wrapper> : std::true_type {}; template constexpr bool is_reference_wrapper_v = is_reference_wrapper::value; template auto INVOKE(T Base::*pmf, Derived&& ref, Args&&... args) noexcept(noexcept((std::forward(ref).*pmf)(std::forward(args)...))) -> std::enable_if_t::value && std::is_base_of>::value, decltype((std::forward(ref).*pmf)(std::forward(args)...))> { return (std::forward(ref).*pmf)(std::forward(args)...); } template auto INVOKE(T Base::*pmf, RefWrap&& ref, Args&&... args) noexcept(noexcept((ref.get().*pmf)(std::forward(args)...))) -> std::enable_if_t::value && is_reference_wrapper_v>, decltype((ref.get().*pmf)(std::forward(args)...))> { return (ref.get().*pmf)(std::forward(args)...); } template auto INVOKE(T Base::*pmf, Pointer&& ptr, Args&&... args) noexcept(noexcept(((*std::forward(ptr)).*pmf)(std::forward(args)...))) -> std::enable_if_t::value && !is_reference_wrapper_v> && !std::is_base_of>::value, decltype(((*std::forward(ptr)).*pmf)(std::forward(args)...))> { return ((*std::forward(ptr)).*pmf)(std::forward(args)...); } template auto INVOKE(T Base::*pmd, Derived&& ref) noexcept(noexcept(std::forward(ref).*pmd)) -> std::enable_if_t::value && std::is_base_of>::value, decltype(std::forward(ref).*pmd)> { return std::forward(ref).*pmd; } template auto INVOKE(T Base::*pmd, RefWrap&& ref) noexcept(noexcept(ref.get().*pmd)) -> std::enable_if_t::value && is_reference_wrapper_v>, decltype(ref.get().*pmd)> { return ref.get().*pmd; } template auto INVOKE(T Base::*pmd, Pointer&& ptr) noexcept(noexcept((*std::forward(ptr)).*pmd)) -> std::enable_if_t::value && !is_reference_wrapper_v> && !std::is_base_of>::value, decltype((*std::forward(ptr)).*pmd)> { return (*std::forward(ptr)).*pmd; } template auto INVOKE(F&& f, Args&&... args) noexcept(noexcept(std::forward(f)(std::forward(args)...))) -> std::enable_if_t>::value, decltype(std::forward(f)(std::forward(args)...))> { return std::forward(f)(std::forward(args)...); } } // namespace detail template< class F, class... ArgTypes > auto invoke(F&& f, ArgTypes&&... args) // exception specification for QoI noexcept(noexcept(detail::INVOKE(std::forward(f), std::forward(args)...))) -> decltype(detail::INVOKE(std::forward(f), std::forward(args)...)) { return detail::INVOKE(std::forward(f), std::forward(args)...); } namespace detail { template constexpr decltype(auto) apply_impl(F &&f, Tuple &&t, std::index_sequence) { return std::invoke(std::forward(f), std::get(std::forward(t))...); } } // namespace detail template constexpr decltype(auto) apply(F &&f, Tuple &&t) { - return detail::apply_impl( - std::forward(f), std::forward(t), - std::make_index_sequence>::value>{}); + return detail::apply_impl + (std::forward(f), std::forward(t), + std::make_index_sequence>::value>{}); } } //namespace std #endif /* UTILITIES_H */ diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh index 6ed7e2f..02873e2 100644 --- a/src/materials/material_base.hh +++ b/src/materials/material_base.hh @@ -1,116 +1,119 @@ /** * file material_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials (constitutive models) * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include "common/common.hh" #include "common/field_map_matrixlike.hh" #include "common/field_collection.hh" #ifndef MATERIAL_BASE_H #define MATERIAL_BASE_H namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class MaterialBase { public: //! typedefs for data handled by this interface //! global field collection for system-wide fields, like stress, strain, etc using GFieldCollection_t = FieldCollection; //! field collection for internal variables, such as eigen-strains, //! plastic strains, damage variables, etc, but also for managing which //! pixels the material is responsible for using MFieldCollection_t = FieldCollection; - using StressMap_t = MatrixFieldMap; - using StrainMap_t = StressMap_t; - using TangentMap_t = MatrixFieldMap; + using StressField_t = TensorField; + using StrainField_t = StressField_t; + using TangentField_t = TensorField; using Ccoord = Ccoord_t; //! Default constructor MaterialBase() = delete; //! Construct by name MaterialBase(std::string name); //! Copy constructor MaterialBase(const MaterialBase &other) = delete; //! Move constructor MaterialBase(MaterialBase &&other) noexcept = delete; //! Destructor virtual ~MaterialBase() noexcept = default; //! Copy assignment operator MaterialBase& operator=(const MaterialBase &other) = delete; //! Move assignment operator MaterialBase& operator=(MaterialBase &&other) noexcept = delete; - //! take responsibility for a pixel identified by its cell coordinates - //! TODO: this won't work. for materials with additional info per pixel (as, e.g. for eigenstrain), we need to pass more parameters, so this needs to be a variadic function. + /** + * take responsibility for a pixel identified by its cell coordinates + * WARNING: this won't work for materials with additional info per pixel + * (as, e.g. for eigenstrain), we need to pass more parameters. Materials + * of this tye need to overload add_pixel + */ void add_pixel(const Ccoord & ccord); //! allocate memory, etc - //! TODO: this won't work. for materials with additional info per pixel (see above TODO), we neet to allocate before we know for sure how many pixels the material is responsible for. virtual void initialize(bool stiffness = false) = 0; //! return the materil's name const std::string & get_name() const; //! for static inheritance stuff constexpr static Dim_t sdim() {return DimS;} constexpr static Dim_t mdim() {return DimM;} //! computes stress - virtual void compute_stress(const StrainMap_t & F, - StressMap_t & P, - Formulation form) = 0; - //! computes stress and tangent stiffness - virtual void compute_stress_stiffness(const StrainMap_t & F, - StressMap_t & P, - TangentMap_t & K, + virtual void compute_stresses(const StrainField_t & F, + StressField_t & P, + Formulation form) = 0; + //! computes stress and tangent moduli + virtual void compute_stresses_tangent(const StrainField_t & F, + StressField_t & P, + TangentField_t & K, Formulation form) = 0; protected: //! members const std::string name; MFieldCollection_t internal_fields; private: }; } // muSpectre #endif /* MATERIAL_BASE_H */ diff --git a/src/materials/material_hyper_elastic1.cc b/src/materials/material_hyper_elastic1.cc index b280702..a4efbd4 100644 --- a/src/materials/material_hyper_elastic1.cc +++ b/src/materials/material_hyper_elastic1.cc @@ -1,66 +1,68 @@ /** * file material_hyper_elastic1.cc * * @author Till Junge * * @date 14 Nov 2017 * * @brief Implementation for materialhyperelastic1 * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "materials/material_hyper_elastic1.hh" #include "common/tensor_algebra.hh" #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template MaterialHyperElastic1::MaterialHyperElastic1(std::string name, Real young, Real poisson) :Parent(name), young{young}, poisson{poisson}, lambda{young*poisson/((1+poisson)*(1-2*poisson))}, mu{young/(2*(1+poisson))}, - C{lambda*Tensors::outer(Tensors::I2(),Tensors::I2())}// + - //2*mu*Tensors::I4S()} + C{lambda*Tensors::outer(Tensors::I2(),Tensors::I2()) + + 2*mu*Tensors::I4S()} {} /* ---------------------------------------------------------------------- */ template + template decltype(auto) - MaterialHyperElastic1::evaluate_stress(const Strain_t & E) { + MaterialHyperElastic1::evaluate_stress(s_t && E) { return E.trace()*lambda * Strain_t::Identity() + 2*mu*E; } /* ---------------------------------------------------------------------- */ template + template decltype(auto) - MaterialHyperElastic1::evaluate_stress_tangent(const Strain_t & E) { - return std::forward_as_tuple(this->evaluate_stress(E), this->C); + MaterialHyperElastic1::evaluate_stress_tangent(s_t && E) { + return std::forward_as_tuple(this->evaluate_stress(std::move(E)), this->C); } template class MaterialHyperElastic1<2, 2>; template class MaterialHyperElastic1<2, 3>; template class MaterialHyperElastic1<3, 3>; } // muSpectre diff --git a/src/materials/material_hyper_elastic1.hh b/src/materials/material_hyper_elastic1.hh index 4e11676..d587aaf 100644 --- a/src/materials/material_hyper_elastic1.hh +++ b/src/materials/material_hyper_elastic1.hh @@ -1,94 +1,101 @@ /** * file material_hyper_elastic1.hh * * @author Till Junge * * @date 13 Nov 2017 * * @brief Implementation for hyperelastic reference material like in de Geus * 2017. This follows the simplest and likely not most efficient * implementation (with exception of the Python law) * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "materials/material_muSpectre_base.hh" #ifndef MATERIAL_HYPER_ELASTIC1_H #define MATERIAL_HYPER_ELASTIC1_H namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class MaterialHyperElastic1: public MaterialMuSpectre, DimS, DimM> { public: using Parent = MaterialMuSpectre; + using GFieldCollection_t = typename Parent::GFieldCollection_t; // declare what type of strain measure your law takes as input constexpr static auto strain_measure{StrainMeasure::GreenLagrange}; // declare what type of stress measure your law yields as output constexpr static auto stress_measure{StressMeasure::PK2}; // declare whether the derivative of stress with respect to strain is uniform constexpr static bool uniform_stiffness = true; // declare the type in which you wish to receive your strain measure - using Strain_t = Eigen::Matrix; - using Stress_t = Strain_t; - using Tangent_t = Eigen::TensorFixedSize - , Eigen::RowMajor>; + using StrainMap_t = MatrixFieldMap; + using StressMap_t = StrainMap_t; + using TangentMap_t = T4MatrixFieldMap; + using Strain_t = typename StrainMap_t::const_reference; + using Stress_t = typename StressMap_t::reference; + using Tangent_t = typename TangentMap_t::reference; + using Stiffness_t = Eigen::TensorFixedSize + >; //! Default constructor MaterialHyperElastic1() = delete; //! Copy constructor MaterialHyperElastic1(const MaterialHyperElastic1 &other) = delete; //! Construct by name, Young's modulus and Poisson's ratio MaterialHyperElastic1(std::string name, Real young, Real poisson); //! Move constructor MaterialHyperElastic1(MaterialHyperElastic1 &&other) noexcept = delete; //! Destructor virtual ~MaterialHyperElastic1() noexcept = default; //! Copy assignment operator MaterialHyperElastic1& operator=(const MaterialHyperElastic1 &other) = delete; //! Move assignment operator MaterialHyperElastic1& operator=(MaterialHyperElastic1 &&other) noexcept = delete; - decltype(auto) evaluate_stress(const Strain_t & E); - decltype(auto) evaluate_stress_tangent(const Strain_t & E); + template + decltype(auto) evaluate_stress(s_t && E); + template + decltype(auto) evaluate_stress_tangent(s_t && E); const Tangent_t & get_stiffness() const; protected: const Real young, poisson, lambda, mu; - const Tangent_t C; + const Stiffness_t C; private: }; } // muSpectre #endif /* MATERIAL_HYPER_ELASTIC1_H */ diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index 1a0af29..e082852 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,358 +1,351 @@ /** * file material_muSpectre_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials written for µSpectre specifically. These * can take full advantage of the configuration-change utilities of * µSpectre. The user can inherit from them to define new constitutive * laws and is merely required to provide the methods for computing the * second Piola-Kirchhoff stress and Tangent. This class uses the * "curiously recurring template parameter" to avoid virtual calls. * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include "stdexcept" #include "materials/material_base.hh" #include "materials/materials_toolbox.hh" #include "common/field_collection.hh" #include "common/field.hh" #ifndef MATERIAL_MUSPECTRE_BASE_H #define MATERIAL_MUSPECTRE_BASE_H namespace muSpectre { //! 'Material' is a CRTP template class MaterialMuSpectre: public MaterialBase { public: enum class NeedTangent {yes, no}; using Parent = MaterialBase; using MFieldCollection_t = typename Parent::MFieldCollection_t; - using StressMap_t = typename Parent::StressMap_t; - using StrainMap_t = typename Parent::StrainMap_t; - using TangentMap_t = typename Parent::TangentMap_t; + using StressField_t = typename Parent::StressField_t; + using StrainField_t = typename Parent::StrainField_t; + using TangentField_t = typename Parent::TangentField_t; //! Default constructor MaterialMuSpectre() = delete; //! Construct by name MaterialMuSpectre(std::string name); //! Copy constructor MaterialMuSpectre(const MaterialMuSpectre &other) = delete; //! Move constructor MaterialMuSpectre(MaterialMuSpectre &&other) noexcept = delete; //! Destructor virtual ~MaterialMuSpectre() noexcept = default; //! Copy assignment operator MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete; //! Move assignment operator MaterialMuSpectre& operator=(MaterialMuSpectre &&other) noexcept = delete; //! computes stress - virtual void compute_stresses(const StrainMap_t & F, - StressMap_t & P, - Formulation form); + virtual void compute_stresses(const StrainField_t & F, + StressField_t & P, + Formulation form) override final; //! computes stress and tangent modulus - virtual void compute_stresses_tangent(const StrainMap_t & F, - StressMap_t & P, - TangentMap_t & K, - Formulation form); - - + virtual void compute_stresses_tangent(const StrainField_t & F, + StressField_t & P, + TangentField_t & K, + Formulation form) override final; protected: //! computes stress with the formulation available at compile time template - inline void compute_stresses_worker(const StrainMap_t & F, - StressMap_t & P); + inline void compute_stresses_worker(const StrainField_t & F, + StressField_t & P); //! computes stress with the formulation available at compile time template - inline void compute_stresses_worker(const StrainMap_t & F, - StressMap_t & P, - TangentMap_t & K); + inline void compute_stresses_worker(const StrainField_t & F, + StressField_t & P, + TangentField_t & K); //! this iterable class is a default for simple laws that just take a strain //! the iterable is just a templated wrapper to provide a range to iterate over //! that does or does not include tangent moduli - template + template class iterable_proxy; private: }; /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: - compute_stresses(const StrainMap_t &F, StressMap_t &P, + compute_stresses(const StrainField_t &F, StressField_t &P, muSpectre::Formulation form) { switch (form) { case Formulation::finite_strain: { this->template compute_stresses_worker(F, P); break; } case Formulation::small_strain: { this->template compute_stresses_worker(F, P); break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: - compute_stresses_tangent(const StrainMap_t & F, StressMap_t & P, - TangentMap_t & K, + compute_stresses_tangent(const StrainField_t & F, StressField_t & P, + TangentField_t & K, muSpectre::Formulation form) { switch (form) { case Formulation::finite_strain: { this->template compute_stresses_worker(F, P, K); break; } case Formulation::small_strain: { this->template compute_stresses_worker(F, P, K); break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre:: - compute_stresses_worker(const StrainMap_t & F, - StressMap_t & P, - TangentMap_t & K){ + compute_stresses_worker(const StrainField_t & F, + StressField_t & P, + TangentField_t & K){ //! in small strain problems, we always want Cauchy stress as a function //! of the infinitesimal strain tensor - using IterableSmallStrain_t = iterable_proxy; - using IterableFiniteStrain_t = iterable_proxy; + using IterableSmallStrain_t = iterable_proxy; + using IterableFiniteStrain_t = iterable_proxy; using iterable = std::conditional_t <(Form==Formulation::small_strain), IterableSmallStrain_t, IterableFiniteStrain_t>; constexpr auto stored_strain_m = get_stored_strain_type(Form); constexpr auto expected_strain_m = get_formulation_strain_type(Form, Material::strain_measure); /* This lambda is executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ auto constitutive_law = [this](auto && F, auto && ... internal_variables) { auto && strain = MatTB::convert_strain(F); // return value contains a tuple of rvalue_refs to both stress and tangent moduli - auto && stress_tgt = this->evaluate_stress_tangent(std::move(strain), internal_variables...); + auto && stress_tgt = static_cast(*this).evaluate_stress_tangent(std::move(strain), internal_variables...); if (Form == Formulation::small_strain) { - return stress_tgt; + return std::move(stress_tgt); } else { auto && stress = std::get<0>(stress_tgt); auto && tangent = std::get<1>(stress_tgt); return MatTB::PK1_stress (F, stress, tangent); } }; - iterable it{this}; + iterable it{*this}; for (auto && tuples: it) { // the iterator yields a pair of tuples. this first tuple contains // references to stress and stiffness in the global arrays, the second // contains references to the deformation gradient and internal variables // (some of them are const). auto && stress_tgt = std::get<0>(tuples); auto && inputs = std::get<1>(tuples); stress_tgt = std::apply(constitutive_law, std::move(inputs)); } } /* ---------------------------------------------------------------------- */ //! this iterator class is a default for simple laws that just take a strain template - template::NeedTangent NeedTgt> + template::NeedTangent NeedTgt> class MaterialMuSpectre::iterable_proxy { public: //! Default constructor iterable_proxy() = delete; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ iterable_proxy(const MaterialMuSpectre & mat); - static constexpr auto strain_measure{StrainMeasure}; - using Strain_t = typename Material::Strain_t::const_reference; - using Stress_t = typename Material::Stress_t::reference; - using Tangent_t = typename Material::Tangent_t::reference; + using Strain_t = typename Material::Strain_t; + using Stress_t = typename Material::Stress_t; + using Tangent_t = typename Material::Tangent_t; //! Copy constructor iterable_proxy(const iterable_proxy &other) = default; //! Move constructor iterable_proxy(iterable_proxy &&other) noexcept = default; //! Destructor virtual ~iterable_proxy() noexcept = default; //! Copy assignment operator iterable_proxy& operator=(const iterable_proxy &other) = default; //! Move assignment operator iterable_proxy& operator=(iterable_proxy &&other) noexcept = default; class iterator { public: using value_type = - std::conditional_t<(NeedTgt == NeedTangent::Yes), - std::tuple, - std::tuple>; + std::conditional_t + <(NeedTgt == NeedTangent::yes), + std::tuple, std::tuple>, + std::tuple, std::tuple>>; using iterator_category = std::forward_iterator_tag; //! Default constructor iterator() = delete; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ iterator(const MaterialMuSpectre & mat, bool begin = true); //! Copy constructor iterator(const iterator &other) = default; //! Move constructor iterator(iterator &&other) noexcept = default; //! Destructor virtual ~iterator() noexcept = default; //! Copy assignment operator iterator& operator=(const iterator &other) = default; //! Move assignment operator iterator& operator=(iterator &&other) noexcept = default; //! pre-increment inline iterator & operator++(); //! dereference inline value_type operator*(); //! inequality inline bool operator!=(const iterator & other) const; protected: const MaterialMuSpectre & material; size_t index; private: }; + iterator begin() {return iterator(*this);} + iterator end() {return iterator(*this, false);} + protected: const MaterialMuSpectre & material; private: }; /* ---------------------------------------------------------------------- */ template - template::NeedTangent Tangent> - MaterialMuSpectre::template iterable_proxy:: + template::NeedTangent Tangent> + MaterialMuSpectre::iterable_proxy:: iterable_proxy(const MaterialMuSpectre & mat) : material{mat} {} /* ---------------------------------------------------------------------- */ template - template::NeedTangent Tangent> - MaterialMuSpectre::iterable_proxy:: + template::NeedTangent Tangent> + MaterialMuSpectre::iterable_proxy:: iterator::iterator(const MaterialMuSpectre & mat, bool begin) : material{mat}, index{begin ? 0:mat.internal_fields.size()} {} /* ---------------------------------------------------------------------- */ //! pre-increment template - template::NeedTangent Tangent> + template::NeedTangent Tangent> typename MaterialMuSpectre:: - template iterable_proxy::iterator & + template iterable_proxy::iterator & MaterialMuSpectre:: - template iterable_proxy::iterator:: + template iterable_proxy::iterator:: operator++() { ++this->index; return *this; } /* ---------------------------------------------------------------------- */ //! dereference template - template::NeedTangent Tangent> + template::NeedTangent Tangent> typename MaterialMuSpectre:: - template iterable_proxy::iterator::value_type + template iterable_proxy::iterator::value_type MaterialMuSpectre:: - template iterable_proxy::iterator:: + template iterable_proxy::iterator:: operator*() { static_assert(Tangent, "Template specialisation failed"); return std::make_tuple (this->material.internal_fields); } } // muSpectre #endif /* MATERIAL_MUSPECTRE_BASE_H */ diff --git a/src/materials/materials_toolbox.hh b/src/materials/materials_toolbox.hh index da2a6a6..dce28e4 100644 --- a/src/materials/materials_toolbox.hh +++ b/src/materials/materials_toolbox.hh @@ -1,306 +1,306 @@ /** * file materials_toolbox.hh * * @author Till Junge * * @date 02 Nov 2017 * * @brief collection of common continuum mechanics tools * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef MATERIALS_TOOLBOX_H #define MATERIALS_TOOLBOX_H #include #include #include #include #include "common/common.hh" #include "common/tensor_algebra.hh" namespace muSpectre { namespace MatTB { class MaterialsToolboxError:public std::runtime_error{ public: explicit MaterialsToolboxError(const std::string& what) :std::runtime_error(what){} explicit MaterialsToolboxError(const char * what) :std::runtime_error(what){} }; /* ---------------------------------------------------------------------- */ /** Structure for functions returning one strain measure as a function of another **/ namespace internal { template struct ConvertStrain { template inline static decltype(auto) compute(Strain_t&& input) { // transparent case, in which no conversion is required: // just a perfect forwarding static_assert ((In == Out), "This particular strain conversion is not implemented"); return std::forward(input); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for getting Green-Lagrange strain from the transformation gradient **/ template <> template decltype(auto) ConvertStrain:: compute(Strain_t && F) { return (F.transpose()*F - Strain_t::PlainObject::Identity()); } } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning one strain measure as a function of //! another template decltype(auto) convert_strain(Strain_t && strain) { return internal::ConvertStrain::compute(std::move(strain)); }; /* ---------------------------------------------------------------------- */ /** Structure for functions returning PK1 stress from other stress measures **/ namespace internal { template struct PK1_stress { - template + template inline static decltype(auto) - compute(Stress_t && /*stress*/, Strain_t && /*strain*/) { + compute(Strain_t && /*strain*/, Stress_t && /*stress*/) { // the following test always fails to generate a compile-time error static_assert((StressM == StressMeasure::Cauchy) && (StressM == StressMeasure::PK1), "The requested Stress conversion is not implemented. " "You either made a programming mistake or need to " "implement it as a specialisation of this function. " "See PK2stress for an example."); } - template + template inline static decltype(auto) - compute(Stress_t && /*stress*/, Tangent_t && /*stiffness*/, - Strain_t && /*strain*/) { + compute(Strain_t && /*strain*/, Stress_t && /*stress*/, + Tangent_t && /*stiffness*/) { // the following test always fails to generate a compile-time error static_assert((StressM == StressMeasure::Cauchy) && (StressM == StressMeasure::PK1), "The requested Stress conversion is not implemented. " "You either made a programming mistake or need to " "implement it as a specialisation of this function. " "See PK2stress for an example."); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for the transparent case, where we already have PK1 stress **/ template struct PK1_stress: public PK1_stress { - template + template inline static decltype(auto) - compute(Stress_t && P, Strain_t && /*dummy*/) { + compute(Strain_t && /*dummy*/, Stress_t && P) { return std::forward(P); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for the transparent case, where we already have PK1 stress *and* stiffness is given with respect to the transformation gradient **/ template <> struct PK1_stress: public PK1_stress { - template + template decltype(auto) - compute(Stress_t && P, Tangent_t && K, Strain_t && /*dummy*/) { + compute(Strain_t && /*dummy*/, Stress_t && P, Tangent_t && K) { return std::forward_as_tuple(P, K); } }; /* ---------------------------------------------------------------------- */ /** * Specialisation for the case where we get material stress (PK2) */ template struct PK1_stress: public PK1_stress { - template + template inline static decltype(auto) - compute(Stress_t && S, Strain_t && F) { + compute(Strain_t && F, Stress_t && S) { return F*S; } }; /* ---------------------------------------------------------------------- */ /** * Specialisation for the case where we get material stress (PK2) derived * with respect to Green-Lagrange strain */ template <> struct PK1_stress: public PK1_stress { - template + template inline static decltype(auto) - compute(Stress_t && S, Tangent_t && C, Strain_t && F) { + compute(Strain_t && F, Stress_t && S, Tangent_t && C) { using T4 = typename Tangent_t::PlainObject; T4 K; K.setZero(); constexpr int dim{Strain_t::ColsAtCompileTime}; for (int i = 0; i < dim; ++i) { for (int m = 0; m < dim; ++m) { for (int n = 0; n < dim; ++n) { K(i,m,i,n) += S(m,n); for (int j = 0; j < dim; ++j) { for (int r = 0; r < dim; ++r) { for (int s = 0; s < dim; ++s) { K(i,m,j,n) += F(i,r)*C(r,m,n,s)*(F(j,s)); } } } } } } return K; } }; } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning an expression for PK2 stress based on template decltype(auto) PK1_stress(Strain_t && strain, Stress_t && stress) { return internal::PK1_stress::compute(std::move(stress), std::move(strain)); }; /* ---------------------------------------------------------------------- */ //! set of functions returning an expression for PK2 stress based on template decltype(auto) PK1_stress(Strain_t && strain, Stress_t && stress, Tangent_t && stiffness) { return internal::PK1_stress::compute(std::move(stress), std::move(strain), std::move(stiffness)); }; /* ---------------------------------------------------------------------- */ /** Structure for functions returning ∂strain/∂F (where F is the transformation * gradient. (e.g. for a material law expressed as PK2 stress S as * function of Green-Lagrange strain E, this would return F^T ⊗̲ I) (⊗_ * refers to Curnier's notation) WARNING: this function assumes that your * stress measure has the two minor symmetries due to stress equilibrium * and the major symmetry due to the existance of an elastic potential * W(E). Do not rely on this function if you are using exotic stress * measures **/ namespace internal { template struct StrainDerivative { template inline static decltype(auto) compute (Grad_t && grad) { // the following test always fails to generate a compile-time error static_assert((StrainM == StrainMeasure::Almansi) && (StrainM == StrainMeasure::Biot), "The requested StrainDerivative calculation is not " "implemented. You either made a programming mistake or " "need to implement it as a specialisation of this " "function. See " "StrainDerivative for an " "example."); } }; template <> template decltype(auto) StrainDerivative:: compute(Grad_t && grad) { } } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning ∂strain/∂F template auto StrainDerivative = [](auto && F) decltype(auto) { // the following test always fails to generate a compile-time error static_assert((StrainM == StrainMeasure::Almansi) && (StrainM == StrainMeasure::Biot), "The requested StrainDerivative calculation is not " "implemented. You either made a programming mistake or " "need to implement it as a specialisation of this " "function. See " "StrainDerivative for an " "example."); }; /* ---------------------------------------------------------------------- */ template<> auto StrainDerivative = [] (auto && F) -> decltype(auto) { constexpr size_t order{2}; constexpr Dim_t dim{F.Dimensions[0]}; return Tensors::outer_under (F.shuffle(std::array{1,0}), Tensors::I2()); }; } // MatTB } // muSpectre #endif /* MATERIALS_TOOLBOX_H */ diff --git a/tests/test_field_collections.cc b/tests/test_field_collections.cc index 09f75b2..d2cf173 100644 --- a/tests/test_field_collections.cc +++ b/tests/test_field_collections.cc @@ -1,548 +1,555 @@ /** * file test_field_collections.cc * * @author Till Junge * * @date 20 Sep 2017 * * @brief Test the FieldCollection classes which provide fast optimized iterators * over run-time typed fields * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include #include #include #include "common/common.hh" #include "common/ccoord_operations.hh" #include "common/test_goodies.hh" #include "tests.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "common/field_map.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(field_collection_tests); //! Test fixture for simple tests on single field in collection template struct FC_fixture: public FieldCollection { FC_fixture() :fc() {} inline static constexpr Dim_t sdim(){return DimS;} inline static constexpr Dim_t mdim(){return DimM;} inline static constexpr bool global(){return Global;} using FC_t = FieldCollection; FC_t fc; }; using test_collections = boost::mpl::list, FC_fixture<2, 3, true>, FC_fixture<3, 3, true>, FC_fixture<2, 2, false>, FC_fixture<2, 3, false>, FC_fixture<3, 3, false>>; BOOST_AUTO_TEST_CASE(simple) { const Dim_t sdim = 2; const Dim_t mdim = 2; using FC_t = FieldCollection; FC_t fc; BOOST_CHECK_EQUAL(FC_t::spatial_dim(), sdim); BOOST_CHECK_EQUAL(fc.get_spatial_dim(), sdim); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(Simple_construction_test, F, test_collections, F) { BOOST_CHECK_EQUAL(F::FC_t::spatial_dim(), F::sdim()); BOOST_CHECK_EQUAL(F::fc.get_spatial_dim(), F::sdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(get_field2_test, F, test_collections, F) { const auto order{2}; using FC_t = typename F::FC_t; using TF_t = TensorField; auto && myfield = make_field("TensorField real 2", F::fc); using TensorMap = TensorFieldMap; using MatrixMap = MatrixFieldMap; using ArrayMap = ArrayFieldMap; TensorMap TFM(myfield); MatrixMap MFM(myfield); ArrayMap AFM(myfield); BOOST_CHECK_EQUAL(TFM.info_string(), "Tensor(d, "+ std::to_string(order) + "_o, " + std::to_string(F::mdim()) + "_d)"); BOOST_CHECK_EQUAL(MFM.info_string(), "Matrix(d, "+ std::to_string(F::mdim()) + "x" + std::to_string(F::mdim()) + ")"); BOOST_CHECK_EQUAL(AFM.info_string(), "Array(d, "+ std::to_string(F::mdim()) + "x" + std::to_string(F::mdim()) + ")"); } constexpr Dim_t order{4}, matrix_order{2}; //! Test fixture for multiple fields in the collection template struct FC_multi_fixture{ FC_multi_fixture() :fc() { //add Real tensor field make_field> ("Tensorfield Real o4", fc); //add Real tensor field make_field> ("Tensorfield Real o2", fc); //add integer scalar field make_field> ("integer Scalar", fc); //add complex matrix field make_field> ("Matrixfield Complex sdim x mdim", fc); } inline static constexpr Dim_t sdim(){return DimS;} inline static constexpr Dim_t mdim(){return DimM;} inline static constexpr bool global(){return Global;} using FC_t = FieldCollection; FC_t fc; }; using mult_collections = boost::mpl::list, FC_multi_fixture<2, 3, true>, FC_multi_fixture<3, 3, true>, FC_multi_fixture<2, 2, false>, FC_multi_fixture<2, 3, false>, FC_multi_fixture<3, 3, false>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(multi_field_test, F, mult_collections, F) { using FC_t = typename F::FC_t; // possible maptypes for Real tensor fields using T_type = Real; using T_TFM1_t = TensorFieldMap; using T_TFM2_t = TensorFieldMap; //! dangerous + using T4_Map_t = T4MatrixFieldMap; + // impossible maptypes for Real tensor fields using T_SFM_t = ScalarFieldMap; using T_MFM_t = MatrixFieldMap; using T_AFM_t = ArrayFieldMap; using T_MFMw1_t = MatrixFieldMap; using T_MFMw2_t = MatrixFieldMap; using T_MFMw3_t = MatrixFieldMap; const std::string T_name{"Tensorfield Real o4"}; const std::string T_name_w{"TensorField Real o4 wrongname"}; BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_NO_THROW(T_TFM1_t(F::fc.at(T_name))); BOOST_CHECK_NO_THROW(T_TFM2_t(F::fc.at(T_name))); + BOOST_CHECK_NO_THROW(T4_Map_t(F::fc.at(T_name))); + BOOST_CHECK_THROW(T4_Map_t(F::fc.at(T_name_w)), std::out_of_range); BOOST_CHECK_THROW(T_MFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_AFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw1_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw3_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name_w)), std::out_of_range); // possible maptypes for integer scalar fields using S_type = Int; using S_SFM_t = ScalarFieldMap; using S_TFM1_t = TensorFieldMap; using S_TFM2_t = TensorFieldMap; using S_MFM_t = MatrixFieldMap; using S_AFM_t = ArrayFieldMap; + using S4_Map_t = T4MatrixFieldMap; // impossible maptypes for integer scalar fields using S_MFMw1_t = MatrixFieldMap; using S_MFMw2_t = MatrixFieldMap; using S_MFMw3_t = MatrixFieldMap; const std::string S_name{"integer Scalar"}; const std::string S_name_w{"integer Scalar wrongname"}; BOOST_CHECK_NO_THROW(S_SFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_TFM1_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_TFM2_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_MFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_AFM_t(F::fc.at(S_name))); + BOOST_CHECK_NO_THROW(S4_Map_t(F::fc.at(S_name))); BOOST_CHECK_THROW(S_MFMw1_t(F::fc.at(S_name)), FieldInterpretationError); + BOOST_CHECK_THROW(T4_Map_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw3_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_SFM_t(F::fc.at(S_name_w)), std::out_of_range); // possible maptypes for complex matrix fields using M_type = Complex; using M_MFM_t = MatrixFieldMap; using M_AFM_t = ArrayFieldMap; // impossible maptypes for complex matrix fields using M_SFM_t = ScalarFieldMap; using M_MFMw1_t = MatrixFieldMap; using M_MFMw2_t = MatrixFieldMap; using M_MFMw3_t = MatrixFieldMap; const std::string M_name{"Matrixfield Complex sdim x mdim"}; const std::string M_name_w{"Matrixfield Complex sdim x mdim wrongname"}; BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_NO_THROW(M_MFM_t(F::fc.at(M_name))); BOOST_CHECK_NO_THROW(M_AFM_t(F::fc.at(M_name))); BOOST_CHECK_THROW(M_MFMw1_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw3_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name_w)), std::out_of_range); } /* ---------------------------------------------------------------------- */ //! Check whether fields can be initialized using mult_collections_t = boost::mpl::list, FC_multi_fixture<2, 3, true>, FC_multi_fixture<3, 3, true>>; using mult_collections_f = boost::mpl::list, FC_multi_fixture<2, 3, false>, FC_multi_fixture<3, 3, false>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_glob, F, mult_collections_t, F) { Ccoord_t size; for (auto && s: size) { s = 3; } BOOST_CHECK_NO_THROW(F::fc.initialise(size)); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca, F, mult_collections_f, F) { testGoodies::RandRange rng; for (int i = 0; i < 7; ++i) { Ccoord_t pixel; for (auto && s: pixel) { s = rng.randval(0, 7); } F::fc.add_pixel(pixel); } BOOST_CHECK_NO_THROW(F::fc.initialise()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca_with_push_back, F, mult_collections_f, F) { constexpr auto mdim = F::mdim(); constexpr int nb_pix = 7; testGoodies::RandRange rng; using ftype = internal::TypedFieldBase; using stype = Eigen::Array; auto & field = reinterpret_cast(F::fc["Tensorfield Real o4"]); field.push_back(stype()); for (int i = 0; i < nb_pix; ++i) { Ccoord_t pixel; for (auto && s: pixel) { s = rng.randval(0, 7); } F::fc.add_pixel(pixel); } BOOST_CHECK_THROW(F::fc.initialise(), FieldCollectionError); for (int i = 0; i < nb_pix-1; ++i) { field.push_back(stype()); } BOOST_CHECK_NO_THROW(F::fc.initialise()); } //! Test fixture for iterators over multiple fields template struct FC_iterator_fixture : public FC_multi_fixture { using Parent = FC_multi_fixture; FC_iterator_fixture() :Parent() { this-> fill(); } template std::enable_if_t fill() { static_assert(Global==isGlobal, "You're breaking my SFINAE plan"); Ccoord_t size; for (auto && s: size) { s = cube_size(); } this->fc.initialise(size); } template std::enable_if_t fill (int dummy = 0) { static_assert(notGlobal != Global, "You're breaking my SFINAE plan"); testGoodies::RandRange rng; this->fc.add_pixel({0,0}); for (int i = 0*dummy; i < sele_size(); ++i) { Ccoord_t pixel; for (auto && s: pixel) { s = rng.randval(0, 7); } this->fc.add_pixel(pixel); } this->fc.initialise(); } constexpr static Dim_t cube_size() {return 3;} constexpr static Dim_t sele_size() {return 7;} }; using iter_collections = boost::mpl::list, FC_iterator_fixture<2, 3, true>, FC_iterator_fixture<3, 3, true>, FC_iterator_fixture<2, 2, false>, FC_iterator_fixture<2, 3, false>, FC_iterator_fixture<3, 3, false>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(iter_field_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using Tensor4Map = TensorFieldMap; Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; F::fc["Tensorfield Real o4"].set_zero(); for (auto && tens:T4map) { BOOST_CHECK_EQUAL(Real(Eigen::Tensor(tens.abs().sum().eval())()), 0); } using Tensor2Map = TensorFieldMap; using MSqMap = MatrixFieldMap; using ASqMap = ArrayFieldMap; Tensor2Map T2map{F::fc["Tensorfield Real o2"]}; MSqMap Mmap{F::fc["Tensorfield Real o2"]}; ASqMap Amap{F::fc["Tensorfield Real o2"]}; auto t2_it = T2map.begin(); auto t2_it_end = T2map.end(); auto m_it = Mmap.begin(); auto a_it = Amap.begin(); for (; t2_it != t2_it_end; ++t2_it, ++m_it, ++a_it) { t2_it->setRandom(); auto && m = *m_it; bool comp = (m == a_it->matrix()); BOOST_CHECK(comp); } using ScalarMap = ScalarFieldMap; ScalarMap s_map{F::fc["integer Scalar"]}; for (Uint i = 0; i < s_map.size(); ++i) { s_map[i] = i; } Uint counter{0}; for (const auto& val: s_map) { BOOST_CHECK_EQUAL(counter++, val); } } using glob_iter_colls = boost::mpl::list, FC_iterator_fixture<2, 3, true>, FC_iterator_fixture<3, 3, true>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(ccoord_indexing_test, F, glob_iter_colls, F) { using FC_t = typename F::Parent::FC_t; using ScalarMap = ScalarFieldMap; ScalarMap s_map{F::fc["integer Scalar"]}; for (Uint i = 0; i < s_map.size(); ++i) { s_map[i] = i; } for (size_t i = 0; i < CcoordOps::get_size(F::fc.get_sizes()); ++i) { Ccoord_t sizes(F::fc.get_sizes()); Ccoord_t ccoords(CcoordOps::get_ccoord(sizes, i)); BOOST_CHECK_EQUAL(CcoordOps::get_index(F::fc.get_sizes(), CcoordOps::get_ccoord(F::fc.get_sizes(), i)), i); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(iterator_methods_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using Tensor4Map = TensorFieldMap; Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; using it_t = typename Tensor4Map::iterator; std::ptrdiff_t diff{3}; // arbitrary, as long as it is smaller than the container size // check constructors auto itstart = T4map.begin(); // standard way of obtaining iterator auto itend = T4map.end(); // ditto it_t it1{T4map}; it_t it2{T4map, false}; it_t it3{T4map, size_t(diff)}; BOOST_CHECK(itstart == itstart); BOOST_CHECK(itstart != itend); BOOST_CHECK_EQUAL(itstart, it1); BOOST_CHECK_EQUAL(itend, it2); // check ostream operator std::stringstream response; response << it3; BOOST_CHECK_EQUAL (response.str(), std::string ("iterator on field 'Tensorfield Real o4', entry ") + std::to_string(diff)); // check copy, move, and assigment constructor (and operator+) it_t itcopy = it1; it_t itmove = std::move(T4map.begin()); it_t it4 = it1+diff; it_t it5(it2); it_t tmp(it2); it_t it6(std::move(tmp)); it_t it7 = it4 -diff; BOOST_CHECK_EQUAL(itcopy, it1); BOOST_CHECK_EQUAL(itmove, it1); BOOST_CHECK_EQUAL(it4, it3); BOOST_CHECK_EQUAL(it5, it2); BOOST_CHECK_EQUAL(it6, it5); BOOST_CHECK_EQUAL(it7, it1); // check increments/decrements BOOST_CHECK_EQUAL(it1++, itcopy); // post-increment BOOST_CHECK_EQUAL(it1, itcopy+1); BOOST_CHECK_EQUAL(--it1, itcopy); // pre -decrement BOOST_CHECK_EQUAL(++it1, itcopy+1); // pre -increment BOOST_CHECK_EQUAL(it1--, itcopy+1); // post-decrement BOOST_CHECK_EQUAL(it1, itcopy); // dereference and member-of-pointer check Eigen::Tensor Tens = *it1; Eigen::Tensor Tens2 = *itstart; Eigen::Tensor check = (Tens==Tens2).all(); BOOST_CHECK_EQUAL(bool(check()), true); BOOST_CHECK_NO_THROW(itstart->setZero()); //check access subscripting auto T3a = *it3; auto T3b = itstart[diff]; BOOST_CHECK(bool(Eigen::Tensor((T3a==T3b).all())())); // div. comparisons BOOST_CHECK_LT(itstart, itend); BOOST_CHECK(!(itend < itstart)); BOOST_CHECK(!(itstart < itstart)); BOOST_CHECK_LE(itstart, itend); BOOST_CHECK_LE(itstart, itstart); BOOST_CHECK(!(itend <= itstart)); BOOST_CHECK_GT(itend, itstart); BOOST_CHECK(!(itend>itend)); BOOST_CHECK(!(itstart>itend)); BOOST_CHECK_GE(itend, itstart); BOOST_CHECK_GE(itend, itend); BOOST_CHECK(!(itstart >= itend)); // check assignment increment/decrement BOOST_CHECK_EQUAL(it1+=diff, it3); BOOST_CHECK_EQUAL(it1-=diff, itstart); // check cell coordinates using Ccoord = Ccoord_t; Ccoord a{itstart.get_ccoord()}; Ccoord b{0}; // Weirdly, boost::has_left_shift::value is false for Ccoords, even though the operator is implemented :( //BOOST_CHECK_EQUAL(a, b); bool check2 = (a==b); BOOST_CHECK(check2); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_tensor_iter_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using Tensor4Map = TensorFieldMap; Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; using T_t = typename Tensor4Map::T_t; Eigen::TensorMap Tens2(T4map[0].data(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim()); for (auto it = T4map.cbegin(); it != T4map.cend(); ++it) { // maps to const tensors can't be initialised with a const pointer this sucks auto&& tens = *it; auto&& ptr = tens.data(); static_assert(std::is_pointer>::value, "should be getting a pointer"); //static_assert(std::is_const>::value, "should be const"); // If Tensor were written well, above static_assert should pass, and the // following check shouldn't. If it get triggered, it means that a newer // version of Eigen now does have const-correct // TensorMap. This means that const-correct field maps // are then also possible for tensors BOOST_CHECK(!std::is_const>::value); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_matrix_iter_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using MatrixMap = MatrixFieldMap; MatrixMap Mmap{F::fc["Matrixfield Complex sdim x mdim"]}; for (auto it = Mmap.cbegin(); it != Mmap.cend(); ++it) { // maps to const tensors can't be initialised with a const pointer this sucks auto&& mat = *it; auto&& ptr = mat.data(); static_assert(std::is_pointer>::value, "should be getting a pointer"); //static_assert(std::is_const>::value, "should be const"); // If Matrix were written well, above static_assert should pass, and the // following check shouldn't. If it get triggered, it means that a newer // version of Eigen now does have const-correct // MatrixMap. This means that const-correct field maps // are then also possible for matrices BOOST_CHECK(!std::is_const>::value); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_scalar_iter_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using ScalarMap = ScalarFieldMap; ScalarMap Smap{F::fc["integer Scalar"]}; for (auto it = Smap.cbegin(); it != Smap.cend(); ++it) { auto&& scal = *it; static_assert(std::is_const>::value, "referred type should be const"); static_assert(std::is_lvalue_reference::value, "Should have returned an lvalue ref"); } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_t4_map.cc b/tests/test_t4_map.cc new file mode 100644 index 0000000..6007092 --- /dev/null +++ b/tests/test_t4_map.cc @@ -0,0 +1,90 @@ +/** + * file test_t4_map.cc + * + * @author Till Junge + * + * @date 20 Nov 2017 + * + * @brief Test the fourth-order map on second-order tensor implementation + * + * @section LICENCE + * + * Copyright (C) 2017 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Emacs; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#include +#include +#include +#include + +#include "common/common.hh" +#include "tests.hh" +#include "common/T4_map_proxy.hh" + + +namespace muSpectre { + + BOOST_AUTO_TEST_SUITE(T4map_tests); + + /** + * Test fixture for construction of T4Map for the time being, symmetry is not + * exploited + */ + template + struct T4_fixture + { + T4_fixture():matrix{}, tensor(matrix.data()){} + + EIGEN_MAKE_ALIGNED_OPERATOR_NEW; + using M4 = Eigen::Matrix; + using T4 = T4Map; + constexpr static Dim_t dim{Dim}; + M4 matrix; + T4 tensor; + }; + + using fix_collection = boost::mpl::list, + T4_fixture>; + + BOOST_FIXTURE_TEST_CASE_TEMPLATE(Simple_construction_test, F, fix_collection, F) { + BOOST_CHECK_EQUAL(F::tensor.cols(), F::dim*F::dim); + } + + BOOST_FIXTURE_TEST_CASE_TEMPLATE(write_access_test, F, fix_collection, F) { + auto & t4 = F::tensor; + constexpr Dim_t dim{F::dim}; + Eigen::TensorFixedSize> t4c; + Eigen::Map t4c_map(t4c.data()); + for (Dim_t i = 0; i < F::dim; ++i) { + for (Dim_t j = 0; j < F::dim; ++j) { + for (Dim_t k = 0; k < F::dim; ++k) { + for (Dim_t l = 0; l < F::dim; ++l) { + t4(i,j,k,l) = 1000*(i+1) + 100*(j+1) + 10*(k+1) + l+1; + t4c(i,j,k,l) = 1000*(i+1) + 100*(j+1) + 10*(k+1) + l+1; + } + } + } + } + for (Dim_t i = 0; i < ipow(dim,4); ++i) { + BOOST_CHECK_EQUAL(F::matrix.data()[i], t4c.data()[i]); + } + } + + BOOST_AUTO_TEST_SUITE_END(); + +} // muSpectre