diff --git a/src/common/geometry.hh b/src/common/geometry.hh index b036eba..eaaeb4f 100644 --- a/src/common/geometry.hh +++ b/src/common/geometry.hh @@ -1,502 +1,502 @@ /** * @file geometry.hh * * @author Till Junge * * @date 18 Apr 2018 * * @brief Geometric calculation helpers * * Copyright © 2018 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 "common/common.hh" #include "common/tensor_algebra.hh" #include "common/eigen_tools.hh" #include #include #include #ifndef SRC_COMMON_GEOMETRY_HH_ #define SRC_COMMON_GEOMETRY_HH_ namespace muSpectre { /** * The rotation matrices depend on the order in which we rotate * around different axes. See [[ * https://en.wikipedia.org/wiki/Euler_angles#Rotation_matrix ]] to * find the matrices */ enum class RotationOrder { Z, XZXEuler, XYXEuler, YXYEuler, YZYEuler, ZYZEuler, ZXZEuler, XZYTaitBryan, XYZTaitBryan, YXZTaitBryan, YZXTaitBryan, ZYXTaitBryan, ZXYTaitBryan }; namespace internal { template struct DefaultOrder { constexpr static RotationOrder value{RotationOrder::ZXYTaitBryan}; }; template <> struct DefaultOrder { constexpr static RotationOrder value{RotationOrder::Z}; }; } // namespace internal template class RotatorBase { public: using RotMat_t = Eigen::Matrix; //! Default constructor RotatorBase() = delete; explicit RotatorBase(const RotMat_t rotation_matrix_input) : rot_mat{rotation_matrix_input} {} //! Copy constructor RotatorBase(const RotatorBase & other) = default; //! Move constructor RotatorBase(RotatorBase && other) = default; //! Destructor virtual ~RotatorBase() = default; //! Copy assignment operator RotatorBase & operator=(const RotatorBase & other) = default; //! Move assignment operator RotatorBase & operator=(RotatorBase && other) = default; /** * Applies the rotation into the frame defined by the rotation * matrix * * @param input is a first-, second-, or fourth-rank tensor * (column vector, square matrix, or T4Matrix, or a Eigen::Map of * either of these, or an expression that evaluates into any of * these) */ template inline decltype(auto) rotate(In_t && input); /** * Applies the rotation back out from the frame defined by the * rotation matrix * * @param input is a first-, second-, or fourth-rank tensor * (column vector, square matrix, or T4Matrix, or a Eigen::Map of * either of these, or an expression that evaluates into any of * these) */ template inline decltype(auto) rotate_back(In_t && input); const RotMat_t & get_rot_mat() const { return this->rot_mat; } void set_rot_mat(const Eigen::Ref & mat_inp) { this->rot_mat = mat_inp; } protected: EIGEN_MAKE_ALIGNED_OPERATOR_NEW; RotMat_t rot_mat; }; template ::value> class RotatorAngle : public RotatorBase { static_assert(((Dim == twoD) and(Order == RotationOrder::Z)) or ((Dim == threeD) and(Order != RotationOrder::Z)), "In 2d, only order 'Z' makes sense. In 3d, it doesn't"); public: using Parent = RotatorBase; using Angles_t = Eigen::Matrix; using RotMat_t = Eigen::Matrix; //! Default constructor RotatorAngle() = delete; explicit RotatorAngle(const Eigen::Ref & angles_inp) : Parent(this->compute_rotation_matrix_angle(angles_inp)), angles{angles_inp} {} //! Copy constructor RotatorAngle(const RotatorAngle & other) = default; //! Move constructor RotatorAngle(RotatorAngle && other) = default; //! Destructor virtual ~RotatorAngle() = default; //! Copy assignment operator RotatorAngle & operator=(const RotatorAngle & other) = default; //! Move assignment operator RotatorAngle & operator=(RotatorAngle && other) = default; protected: inline RotMat_t compute_rotation_matrix_angle(Angles_t angles); inline RotMat_t compute_this_rotation_matrix_angle() { return compute_rotation_matrix(this->angles); } EIGEN_MAKE_ALIGNED_OPERATOR_NEW; Angles_t angles; private: }; /* ---------------------------------------------------------------------- */ namespace internal { template struct RotationMatrixComputerAngle {}; template struct RotationMatrixComputerAngle { constexpr static Dim_t Dim{twoD}; using RotMat_t = typename RotatorAngle::RotMat_t; using Angles_t = typename RotatorAngle::Angles_t; inline static decltype(auto) compute(const Eigen::Ref & angles) { static_assert(Order == RotationOrder::Z, "Two-d rotations can only be around the z axis"); return RotMat_t(Eigen::Rotation2Dd(angles(0))); } }; template struct RotationMatrixComputerAngle { constexpr static Dim_t Dim{threeD}; using RotMat_t = typename RotatorAngle::RotMat_t; using Angles_t = typename RotatorAngle::Angles_t; inline static decltype(auto) compute(const Eigen::Ref & angles) { static_assert(Order != RotationOrder::Z, "three-d rotations cannot only be around the z axis"); switch (Order) { case RotationOrder::ZXZEuler: { return RotMat_t( (Eigen::AngleAxisd(angles(0), Eigen::Vector3d::UnitZ()) * Eigen::AngleAxisd(angles(1), Eigen::Vector3d::UnitX()) * Eigen::AngleAxisd(angles(2), Eigen::Vector3d::UnitZ()))); break; } case RotationOrder::ZXYTaitBryan: { return RotMat_t( (Eigen::AngleAxisd(angles(0), Eigen::Vector3d::UnitZ()) * Eigen::AngleAxisd(angles(1), Eigen::Vector3d::UnitX()) * Eigen::AngleAxisd(angles(2), Eigen::Vector3d::UnitY()))); } default: { throw std::runtime_error("not yet implemented."); } } } }; } // namespace internal /* ---------------------------------------------------------------------- */ template auto RotatorAngle::compute_rotation_matrix_angle(Angles_t angles) -> RotMat_t { return internal::RotationMatrixComputerAngle::compute(angles); } /* ---------------------------------------------------------------------- */ /** * this class is used to make the vector a aligned to the vec b by means of a rotation system, the input for the constructor is the vector itself and the functions rotate and rotate back would be available as they exist in the parent class - (RotatorBase) nad can be used in order to do the functionality of the class + (RotatorBase) and can be used in order to do the functionality of the class */ template class RotatorTwoVec : public RotatorBase { public: using Parent = RotatorBase; using Vec_t = Eigen::Matrix; using RotMat_t = Eigen::Matrix; //! Default constructor RotatorTwoVec() = delete; RotatorTwoVec(Vec_t vec_a_inp, Vec_t vec_b_inp) : Parent(this->compute_rotation_matrix_TwoVec(vec_a_inp, vec_b_inp)), vec_ref{vec_a_inp}, vec_des{vec_b_inp} {} //! Copy constructor RotatorTwoVec(const RotatorTwoVec & other) = default; //! Move constructor RotatorTwoVec(RotatorTwoVec && other) = default; //! Destructor virtual ~RotatorTwoVec() = default; //! Copy assignment operator RotatorTwoVec & operator=(const RotatorTwoVec & other) = default; //! Move assignment operator RotatorTwoVec & operator=(RotatorTwoVec && other) = default; protected: inline RotMat_t compute_rotation_matrix_TwoVec(Vec_t vec_ref, Vec_t vec_des); inline RotMat_t compute_this_rotation_matrix_TwoVec() { return compute_rotation_matrix(this->vec_ref, this->vec_des); } EIGEN_MAKE_ALIGNED_OPERATOR_NEW; Vec_t vec_ref, vec_des; private: }; /* ---------------------------------------------------------------------- */ namespace internal { template struct RotationMatrixComputerTwoVec {}; template <> struct RotationMatrixComputerTwoVec { constexpr static Dim_t Dim{twoD}; using RotMat_t = typename RotatorTwoVec::RotMat_t; using Vec_t = typename RotatorTwoVec::Vec_t; inline static RotMat_t compute(Vec_t vec_ref, Vec_t vec_des) { Real v_ref_norm = sqrt(vec_ref(0) * vec_ref(0) + vec_ref(1) * vec_ref(1)); Real v_des_norm = sqrt(vec_des(0) * vec_des(0) + vec_des(1) * vec_des(1)); RotMat_t ret_mat; ret_mat(0, 0) = ret_mat(1, 1) = (((vec_ref(0) / v_ref_norm) * (vec_des(0) / v_des_norm)) + ((vec_des(1) / v_des_norm) * (vec_ref(1) / v_ref_norm))); ret_mat(1, 0) = (((vec_ref(0) / v_ref_norm) * (vec_des(1) / v_des_norm)) - ((vec_des(0) / v_des_norm) * (vec_ref(1) / v_ref_norm))); ret_mat(0, 1) = -ret_mat(1, 0); return ret_mat; } }; template <> struct RotationMatrixComputerTwoVec { constexpr static Dim_t Dim{threeD}; using RotMat_t = typename RotatorTwoVec::RotMat_t; using Vec_t = typename RotatorTwoVec::Vec_t; inline static RotMat_t compute(Vec_t vec_ref, Vec_t vec_des) { return Eigen::Quaternion::FromTwoVectors(vec_ref, vec_des) .normalized() .toRotationMatrix(); } }; } // namespace internal /* ---------------------------------------------------------------------- */ template auto RotatorTwoVec::compute_rotation_matrix_TwoVec(Vec_t vec_ref, Vec_t vec_des) -> RotMat_t { return internal::RotationMatrixComputerTwoVec::compute(vec_ref, vec_des); } /* ---------------------------------------------------------------------- */ /** * this class is used to make a vector aligned to x-axis of the coordinate system, the input for the constructor is the vector itself and the functions rotate and rotate back would be available as they exist in the parent class (RotatorBase) nad can be used in order to do the functionality of the class */ template class RotatorNormal : public RotatorBase { public: using Parent = RotatorBase; using Vec_t = Eigen::Matrix; using RotMat_t = Eigen::Matrix; //! Default constructor RotatorNormal() = delete; explicit RotatorNormal(Vec_t vec) : Parent(this->compute_rotation_matrix_normal(vec)), vec{vec} {} //! Copy constructor RotatorNormal(const RotatorNormal & other) = default; //! Move constructor RotatorNormal(RotatorNormal && other) = default; //! Destructor virtual ~RotatorNormal() = default; //! Copy assignment operator RotatorNormal & operator=(const RotatorNormal & other) = default; //! Move assignment operator RotatorNormal & operator=(RotatorNormal && other) = default; protected: inline RotMat_t compute_rotation_matrix_normal(Vec_t vec); inline RotMat_t compute_this_rotation_matrix_normal() { return compute_rotation_matrix_normal(this->vec); } EIGEN_MAKE_ALIGNED_OPERATOR_NEW; Vec_t vec; private: }; /* ---------------------------------------------------------------------- */ namespace internal { template struct RotationMatrixComputerNormal {}; template <> struct RotationMatrixComputerNormal { constexpr static Dim_t Dim{twoD}; using RotMat_t = typename RotatorTwoVec::RotMat_t; using Vec_t = typename RotatorTwoVec::Vec_t; inline static RotMat_t compute(Vec_t vec) { Real vec_norm = sqrt(vec(0) * vec(0) + vec(1) * vec(1)); Vec_t x; x << 1.0, 0.0; RotMat_t ret_mat; ret_mat(0, 0) = ret_mat(1, 1) = ((vec(0) / vec_norm) * x(0)); ret_mat(1, 0) = -(-(vec(1) / vec_norm) * x(0)); ret_mat(0, 1) = -ret_mat(1, 0); return ret_mat; } }; template <> struct RotationMatrixComputerNormal { constexpr static Dim_t Dim{threeD}; using RotMat_t = typename RotatorTwoVec::RotMat_t; using Vec_t = typename RotatorTwoVec::Vec_t; inline static RotMat_t compute(Vec_t vec) { Real eps = 0.1; Vec_t vec1 = vec / vec.norm(); Vec_t x(Vec_t::UnitX()); Vec_t y(Vec_t::UnitY()); Vec_t n_x = vec1.cross(x); Vec_t vec2 = ((n_x.norm() > eps) * n_x + (1 - (n_x.norm() > eps)) * (vec1.cross(y))); Vec_t vec3 = vec1.cross(vec2); RotMat_t ret_mat; ret_mat << vec1(0), vec2(0) / vec2.norm(), vec3(0) / vec3.norm(), vec1(1), vec2(1) / vec2.norm(), vec3(1) / vec3.norm(), vec1(2), vec2(2) / vec2.norm(), vec3(2) / vec3.norm(); return ret_mat; } }; } // namespace internal /* ---------------------------------------------------------------------- */ template auto RotatorNormal::compute_rotation_matrix_normal(Vec_t vec) -> RotMat_t { return internal::RotationMatrixComputerNormal::compute(vec); } /* ---------------------------------------------------------------------- */ namespace internal { template struct RotationHelper {}; /* ---------------------------------------------------------------------- */ template <> struct RotationHelper { template inline static decltype(auto) rotate(In_t && input, Rot_t && R) { return R * input; } }; /* ---------------------------------------------------------------------- */ template <> struct RotationHelper { template inline static decltype(auto) rotate(In_t && input, Rot_t && R) { return R * input * R.transpose(); } }; /* ---------------------------------------------------------------------- */ template <> struct RotationHelper { template inline static decltype(auto) rotate(In_t && input, Rot_t && R) { constexpr Dim_t Dim{EigenCheck::tensor_dim::value}; auto && rotator_forward{ Matrices::outer_under(R.transpose(), R.transpose())}; auto && rotator_back = Matrices::outer_under(R, R); // unclear behaviour. When I return this value as an // expression, clange segfaults or returns an uninitialised // tensor return T4Mat(rotator_back * input * rotator_forward); } }; } // namespace internal /* ---------------------------------------------------------------------- */ template template auto RotatorBase::rotate(In_t && input) -> decltype(auto) { constexpr Dim_t tensor_rank{EigenCheck::tensor_rank::value}; return internal::RotationHelper::rotate( std::forward(input), this->rot_mat); } /* ---------------------------------------------------------------------- */ template template auto RotatorBase::rotate_back(In_t && input) -> decltype(auto) { constexpr Dim_t tensor_rank{EigenCheck::tensor_rank::value}; return internal::RotationHelper::rotate( std::forward(input), this->rot_mat.transpose()); } } // namespace muSpectre #endif // SRC_COMMON_GEOMETRY_HH_ diff --git a/src/materials/laminate_homogenisation.cc b/src/materials/laminate_homogenisation.cc index 11754d2..a8db975 100644 --- a/src/materials/laminate_homogenisation.cc +++ b/src/materials/laminate_homogenisation.cc @@ -1,626 +1,627 @@ /** * @file lamminate_hemogenization.hh * * @author Ali Falsafi * * @date 28 Sep 2018 * * @brief * * Copyright © 2017 Till Junge, Ali Falsafi * * µ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 "common/geometry.hh" #include "common/common.hh" #include "laminate_homogenisation.hh" #include "material_anisotropic.hh" #include "materials_toolbox.hh" #include namespace muSpectre { /*---------------------------------------------------------------------- */ template <> constexpr auto LamHomogen::get_serial_indexes() -> std::array { constexpr Dim_t Dim{3}; return std::array{0, 1, 2, 3, 6}; } template <> constexpr auto LamHomogen::get_serial_indexes() -> std::array { constexpr Dim_t Dim{2}; return std::array{0, 1, 2}; } /* ---------------------------------------------------------------------- */ template <> constexpr auto LamHomogen::get_parallel_indexes() -> std::array { constexpr Dim_t Dim{3}; return std::array{4, 5, 7, 8}; } template <> constexpr auto LamHomogen::get_parallel_indexes() -> std::array { constexpr Dim_t Dim{2}; return std::array{3}; } /* ---------------------------------------------------------------------- */ template auto LamHomogen::get_serial_strain(const Eigen::Ref & Fs) -> Serial_strain_t { Serial_strain_t ret_F; auto serial_indexes{get_serial_indexes()}; for (int i{0}; i < ret_F.size(); ++i) { ret_F(i) = Fs(serial_indexes[i]); } return ret_F; } /* ---------------------------------------------------------------------- */ template auto LamHomogen::get_parallel_strain(const Eigen::Ref & Fs) -> Parallel_strain_t { Parallel_strain_t ret_F; auto parallel_indexes(get_parallel_indexes()); for (int i{0}; i < ret_F.size(); ++i) { ret_F(i) = Fs(parallel_indexes[i]); } return ret_F; } /* ---------------------------------------------------------------------- */ template auto LamHomogen::get_equation_strain_from_serial( const Eigen::Ref & Fs) -> Equation_strain_t { Equation_strain_t ret_F; for (int i{0}; i < ret_F.size(); ++i) { ret_F(i) = Fs(i); } return ret_F; } template auto LamHomogen::get_serial_strain_from_equation( const Eigen::Ref & Fe) -> Serial_strain_t { Serial_strain_t ret_F; for (unsigned int j = 0; j < Dim; ++j) { ret_F(j) = Fe(j); } for (unsigned int j = 1; j < Dim; ++j) { ret_F(j + (Dim - 1)) = Fe(j); } return ret_F; } /* ---------------------------------------------------------------------- */ template auto LamHomogen::get_total_strain( const Eigen::Ref & F_seri, const Eigen::Ref & F_para) -> Strain_t { Strain_t ret_F; auto serial_indexes{get_serial_indexes()}; auto parallel_indexes{get_parallel_indexes()}; for (unsigned int i{0}; i < serial_indexes.size(); ++i) { ret_F(serial_indexes[i]) = F_seri(i); } for (unsigned int i{0}; i < parallel_indexes.size(); ++i) { ret_F(parallel_indexes[i]) = F_para(i); } return ret_F; } /* ---------------------------------------------------------------------- */ template auto LamHomogen::get_serial_stiffness(const Eigen::Ref & K) -> Serial_stiffness_t { Serial_stiffness_t ret_K{}; auto serial_indexes{get_serial_indexes()}; for (int i{0}; i < (2 * Dim - 1); ++i) { for (int j{0}; j < (2 * Dim - 1); ++j) { ret_K(i, j) = K(serial_indexes[i], serial_indexes[j]); } } return ret_K; } /* ---------------------------------------------------------------------- */ template auto LamHomogen::get_parallel_stiffness(const Eigen::Ref & K) -> Parallel_stiffness_t { Parallel_stiffness_t ret_K{}; auto parallel_indexes{get_parallel_indexes()}; for (int i{0}; i < (Dim - 1) * (Dim - 1); ++i) { for (int j{0}; j < (Dim - 1) * (Dim - 1); ++j) { ret_K(i, j) = K(parallel_indexes[i], parallel_indexes[j]); } } return ret_K; } /* ---------------------------------------------------------------------- */ template auto LamHomogen::get_total_stiffness( const Eigen::Ref & K_seri, const Eigen::Ref & K_para) -> Stiffness_t { Stiffness_t ret_K{}; auto serial_indexes{get_serial_indexes()}; auto parallel_indexes{get_parallel_indexes()}; for (int i{0}; i < (Dim - 1) * (Dim - 1); ++i) { for (int j{0}; j < (Dim - 1) * (Dim - 1); ++j) { ret_K(parallel_indexes[i], parallel_indexes[j]) = K_para(i, j); } } for (int i{0}; i < (Dim - 1) * (Dim - 1); ++i) { for (int j{0}; j < (Dim - 1) * (Dim - 1); ++j) { ret_K(serial_indexes[i], serial_indexes[j]) = K_seri(i, j); } } return ret_K; } /* ---------------------------------------------------------------------- */ // TODO(falsafi): complete this function and its inverse template auto LamHomogen::get_equation_stiffness_from_serial( const Eigen::Ref & Ks) -> Equation_stiffness_t { Equation_stiffness_t ret_K{}; for (Dim_t i = 0; i < Dim; ++i) { for (Dim_t j = 0; j < Dim; ++j) { if (i == 0 and j == 0) { ret_K(i, j) = 1 * Ks(i, j); } else { ret_K(i, j) = 2 * Ks(i, j); } } } return ret_K; } // namespace muSpectre /* ---------------------------------------------------------------------- */ // It does not to be especialised because all functions used init are // especialised. this function computes the strain in the second // layer of laminate material from strain in the first layer and total strain template auto LamHomogen::linear_eqs_seri(Real ratio, const Eigen::Ref & F_0, const Eigen::Ref & F_1_seri) -> Serial_strain_t { Serial_strain_t F_2_seri; Serial_strain_t F_0_seri = get_serial_strain(F_0); F_2_seri = (F_0_seri - ratio * F_1_seri) / (1 - ratio); return F_2_seri; } /* ---------------------------------------------------------------------- */ // this function recieves the strain in general corrdinate and returns ths // delta_stress and the jacobian in rotated coordinates template auto LamHomogen::delta_stress_eval(const Function_t & mat_1_stress_eval, const Function_t & mat_2_stress_eval, const Eigen::Ref & F_1, const Eigen::Ref & F_2, RotatorNormal rotator, Real ratio) -> std::tuple { auto P_K_1 = mat_1_stress_eval(F_1); auto P_K_2 = mat_2_stress_eval(F_2); StressMat_t del_P_mat; del_P_mat = std::get<0>(P_K_1) - std::get<0>(P_K_2); del_P_mat = rotator.rotate_back(del_P_mat).eval(); auto del_P_seri = get_serial_strain(Eigen::Map(del_P_mat.data())); // P₁(F₁) - P₂(F₂(F₁)) => jacobian of the second term is : ∂P₂/∂F₂ * ∂F₂/∂F₁ Stiffness_t del_K; // std::cout << " K_1 :" << std::endl << std::get<1>(P_K_1) << std::endl; // std::cout << " K_2 :" << std::endl << std::get<1>(P_K_2) << std::endl; del_K = (std::get<1>(P_K_1) - (-ratio / (1 - ratio)) * std::get<1>(P_K_2)); // std::cout << " del_K :" << std::endl << del_K << std::endl; del_K = rotator.rotate_back(del_K).eval(); // std::cout << " del_K :" << std::endl << del_K << std::endl; auto del_K_seri = get_serial_stiffness(del_K); return std::make_tuple(del_P_seri, del_K_seri); } /* ---------------------------------------------------------------------- */ // linear equation are solved in rot coordiantes to obtain F_2 // form known F_1 and then these two are passed to calculate // stress in the general corrdinates template auto LamHomogen::delta_stress_eval_F_1( const Function_t & mat_1_stress_eval, const Function_t & mat_2_stress_eval, Eigen::Ref F_0, Eigen::Ref F_1_rot, RotatorNormal rotator, Real ratio) -> std::tuple { StrainMat_t F_0_rot = rotator.rotate_back(Eigen::Map(F_0.data())); Parallel_strain_t F_0_par_rot = get_parallel_strain(Eigen::Map(F_0_rot.data())); StrainMat_t F_1 = rotator.rotate(Eigen::Map(F_1_rot.data())); Serial_strain_t F_1_seri_rot = get_serial_strain(Eigen::Map(F_1_rot.data())); Serial_strain_t F_2_seri_rot = linear_eqs_seri( ratio, Eigen::Map(F_0_rot.data()), F_1_seri_rot); Strain_t F_2_rot = get_total_strain(F_2_seri_rot, F_0_par_rot); StrainMat_t F_2 = rotator.rotate(Eigen::Map(F_2_rot.data())); return delta_stress_eval(mat_1_stress_eval, mat_2_stress_eval, Eigen::Map(F_1.data()), Eigen::Map(F_2.data()), rotator, ratio); } /* ---------------------------------------------------------------------- */ template inline auto LamHomogen::del_energy_eval(const Real del_F_norm, const Real delta_P_norm) -> Real { return del_F_norm * delta_P_norm; } /* ---------------------------------------------------------------------- */ template <> auto LamHomogen::lam_K_combine(const Eigen::Ref & K_1, const Eigen::Ref & K_2, const Real ratio) -> Stiffness_t { using Mat_A_t = Eigen::Matrix; using Vec_A_t = Eigen::Matrix; using Vec_AT_t = Eigen::Matrix; std::array cf = {1.0, sqrt(2.0), 2.0}; auto get_A11 = [cf](const Eigen::Ref & K) { Mat_A_t A11 = Mat_A_t::Zero(); A11 << cf[0] * get(K, 0, 0, 0, 0), cf[1] * get(K, 0, 0, 0, 1), cf[1] * get(K, 0, 0, 0, 1), cf[2] * get(K, 0, 1, 0, 1); return A11; }; Mat_A_t A11c; A11c << cf[0], cf[1], cf[1], cf[2]; auto get_A12 = [cf](const Eigen::Ref & K) { Vec_A_t A12 = Vec_A_t::Zero(); A12 << cf[0] * get(K, 0, 0, 1, 1), cf[1] * get(K, 1, 1, 0, 1); return A12; }; Vec_A_t A12c = Vec_A_t::Zero(); A12c << cf[0], cf[1]; Mat_A_t A11_1{get_A11(K_1)}; Mat_A_t A11_2{get_A11(K_2)}; Vec_A_t A12_1 = {get_A12(K_1)}; Vec_AT_t A21_1 = A12_1.transpose(); Vec_A_t A12_2 = {get_A12(K_2)}; Vec_AT_t A21_2 = A12_2.transpose(); Real A22_1 = get(K_1, 1, 1, 1, 1); Real A22_2 = get(K_2, 1, 1, 1, 1); auto get_inverse_average = [&ratio](const Eigen::Ref & matrix_1, const Eigen::Ref & matrix_2) { return ((ratio * matrix_1.inverse() + (1 - ratio) * matrix_2.inverse()) .inverse()); }; auto get_average = [&ratio](Real A_1, Real A_2) { return ratio * A_1 + (1 - ratio) * A_2; }; auto get_average_vec = [&ratio](Vec_A_t A_1, Vec_A_t A_2) { return ratio * A_1 + (1 - ratio) * A_2; }; auto get_average_vecT = [&ratio](Vec_AT_t A_1, Vec_AT_t A_2) { return ratio * A_1 + (1 - ratio) * A_2; }; // calculating average of A matrices of the materials Mat_A_t A11 = get_inverse_average(A11_1, A11_2); Vec_A_t A12 = A11 * get_average_vec(A11_1.inverse() * A12_1, A11_2.inverse() * A12_2); Real A22 = get_average(A22_1 - A21_1 * A11_1.inverse() * A12_1, A22_2 - A21_2 * A11_2.inverse() * A12_2) + get_average_vecT(A21_1 * A11_1.inverse(), A21_2 * A11_2.inverse()) * A11 * get_average_vec(A11_1.inverse() * A12_1, A11_2.inverse() * A12_2); std::vector c_maker_inp = { A11(0, 0) / A11c(0, 0), A12(0, 0) / A12c(0, 0), A11(0, 1) / A11c(0, 1), A22, A12(1, 0) / A12c(1, 0), A11(1, 1) / A11c(1, 1)}; // now the resultant stiffness is calculated fro 6 elements obtained from // A matrices averaging routine: Stiffness_t ret_K = Stiffness_t::Zero(); ret_K = MaterialAnisotropic::c_maker(c_maker_inp); return ret_K; } /*------------------------------------------------------------------*/ template <> auto LamHomogen::lam_K_combine(const Eigen::Ref & K_1, const Eigen::Ref & K_2, const Real ratio) -> Stiffness_t { // the combination method is obtained form P. 163 of // "Theory of Composites" // Author : Milton_G_W // constructing "A" matrices( A11, A12, A21, A22) according to the // procedure from the book: // this type of matrix will be used in calculating the combinatio of the // Stiffness matrixes using Mat_A_t = Eigen::Matrix; // these coeffs are used in constructing A matrices from k and vice versa std::array cf = {1.0, sqrt(2.0), 2.0}; auto get_A11 = [cf](const Eigen::Ref & K) { Mat_A_t A11 = Mat_A_t::Zero(); A11 << cf[0] * get(K, 0, 0, 0, 0), cf[1] * get(K, 0, 0, 0, 2), cf[1] * get(K, 0, 0, 0, 1), cf[1] * get(K, 0, 0, 0, 2), cf[2] * get(K, 0, 2, 0, 2), cf[2] * get(K, 0, 2, 0, 1), cf[1] * get(K, 0, 0, 0, 1), cf[2] * get(K, 0, 2, 0, 1), cf[2] * get(K, 0, 1, 0, 1); return A11; }; Mat_A_t A11c; A11c << cf[0], cf[1], cf[1], cf[1], cf[2], cf[2], cf[1], cf[2], cf[2]; auto get_A12 = [cf](const Eigen::Ref & K) { Mat_A_t A12 = Mat_A_t::Zero(); A12 << cf[0] * get(K, 0, 0, 1, 1), cf[0] * get(K, 0, 0, 2, 2), cf[1] * get(K, 0, 0, 1, 2), cf[1] * get(K, 1, 1, 0, 2), cf[1] * get(K, 2, 2, 0, 2), cf[2] * get(K, 1, 2, 0, 2), cf[1] * get(K, 1, 1, 0, 1), cf[1] * get(K, 2, 2, 0, 1), cf[2] * get(K, 1, 2, 0, 1); return A12; }; Mat_A_t A12c; A12c << cf[0], cf[0], cf[1], cf[1], cf[1], cf[2], cf[1], cf[1], cf[2]; auto get_A22 = [cf](const Eigen::Ref & K) { Mat_A_t A22 = Mat_A_t::Zero(); A22 << cf[0] * get(K, 2, 2, 2, 2), cf[0] * get(K, 2, 2, 1, 1), cf[1] * get(K, 2, 2, 1, 2), cf[0] * get(K, 2, 2, 1, 1), cf[0] * get(K, 1, 1, 1, 1), cf[1] * get(K, 1, 1, 1, 2), cf[1] * get(K, 2, 2, 1, 2), cf[1] * get(K, 1, 1, 1, 2), cf[2] * get(K, 2, 1, 2, 1); return A22; }; Mat_A_t A22c; A22c << cf[0], cf[0], cf[1], cf[0], cf[0], cf[1], cf[1], cf[1], cf[2]; Mat_A_t A11_1{get_A11(K_1)}; Mat_A_t A11_2{get_A11(K_2)}; Mat_A_t A12_1{get_A12(K_1)}; Mat_A_t A21_1{A12_1.transpose()}; Mat_A_t A12_2{get_A12(K_2)}; Mat_A_t A21_2{A12_2.transpose()}; Mat_A_t A22_1{get_A22(K_1)}; Mat_A_t A22_2{get_A22(K_2)}; // these two functions are routines to compute average of "A" matrices auto get_inverse_average = [&ratio](const Eigen::Ref & matrix_1, const Eigen::Ref & matrix_2) { return ((ratio * matrix_1.inverse() + (1 - ratio) * matrix_2.inverse()) .inverse()); }; auto get_average = [&ratio](const Mat_A_t & matrix_1, const Mat_A_t & matrix_2) { return (ratio * matrix_1 + (1 - ratio) * matrix_2); }; // calculating average of A matrices of the materials Mat_A_t A11 = get_inverse_average(A11_1, A11_2); Mat_A_t A12 = A11 * get_average(A11_1.inverse() * A12_1, A11_2.inverse() * A12_2); Mat_A_t A22 = get_average(A22_1 - A21_1 * A11_1.inverse() * A12_1, A22_2 - A21_2 * A11_2.inverse() * A12_2) + get_average(A21_1 * A11_1.inverse(), A21_2 * A11_2.inverse()) * A11 * get_average(A11_1.inverse() * A12_1, A11_2.inverse() * A12_2); std::vector c_maker_inp = { A11(0, 0) / A11c(0, 0), A12(0, 0) / A12c(0, 0), A12(0, 1) / A12c(0, 1), A12(0, 2) / A12c(0, 2), A11(1, 0) / A11c(1, 0), A11(0, 2) / A11c(0, 2), A22(1, 1) / A22c(1, 1), A22(1, 0) / A22c(1, 0), A22(2, 1) / A22c(2, 1), A12(1, 0) / A12c(1, 0), A12(2, 0) / A12c(2, 0), A22(0, 0) / A22c(0, 0), A22(2, 0) / A22c(2, 0), A12(1, 1) / A12c(1, 1), A12(2, 1) / A12c(2, 1), A22(2, 2) / A22c(2, 2), A12(1, 2) / A12c(1, 2), A12(2, 2) / A12c(2, 2), A11(1, 1) / A11c(1, 1), A11(2, 1) / A11c(2, 1), A11(2, 2) / A11c(2, 2)}; // now the resultant stiffness is calculated fro 21 elements obtained from // A matrices averaging routine : Stiffness_t ret_K = Stiffness_t::Zero(); ret_K = MaterialAnisotropic::c_maker(c_maker_inp); return ret_K; } /* ---------------------------------------------------------------------- */ template auto LamHomogen::lam_P_combine(const Eigen::Ref & P_1, const Eigen::Ref & P_2, const Real ratio) -> StressMat_t { /*Serial_stress_t P_1_seri = get_serial_strain(Eigen::Map(P_1.data())); Parallel_stress_t P_1_para = get_parallel_strain(Eigen::Map(P_1.data())); Parallel_stress_t P_2_para = get_parallel_strain(Eigen::Map(P_2.data())); Parallel_stress_t P_para = ratio * P_1_para + (1 - ratio) * P_2_para; auto ret_1_P = get_total_strain(P_1_seri, P_para);*/ auto ret_P = ratio * P_1 + (1 - ratio) * P_2; return ret_P; // Eigen::Map(ret_P.data()); } /* ---------------------------------------------------------------------- */ template auto LamHomogen::laminate_solver(Eigen::Ref F_coord, Function_t mat_1_stress_eval, Function_t mat_2_stress_eval, Real ratio, const Eigen::Ref & normal_vec, Real tol, Dim_t max_iter) -> std::tuple { /* * here we rotate the strain such that the laminate intersection normal * would align with the x-axis F_lam is the total strain in the new * coordinates. */ Real del_energy; Dim_t iter{0}; Vec_t my_normal_vec{normal_vec}; RotatorNormal rotator(normal_vec); Strain_t F_1{F_coord}, F_2{F_coord}; StrainMat_t F_mat{rotator.rotate_back(Eigen::Map(F_1.data()))}; auto F_1_mat{F_mat}; auto F_2_mat{F_mat}; Strain_t F_1_rot = Eigen::Map(F_mat.data()); Strain_t F_0_rot = F_1_rot; Parallel_strain_t F_para{get_parallel_strain(F_1_rot)}; Serial_strain_t F_1_new_seri{get_serial_strain(F_1_rot)}; Serial_strain_t F_1_seri{F_1_new_seri}; Equation_strain_t F_1_eq = get_equation_strain_from_serial(F_1_seri); Equation_strain_t F_1_new_eq = get_equation_strain_from_serial(F_1_seri); std::tuple delta_P_K_seri; Serial_stress_t delta_P_seri{}; Serial_stiffness_t delta_K_seri{}; Equation_stress_t delta_P_eq{}; Equation_stiffness_t delta_K_eq{}; std::tuple P_K_1, P_K_2; StressMat_t P_1{}, P_2{}, ret_P{}; Stiffness_t K_1{}, K_2{}, ret_K{}; /* (P_1 - P_2) (as functionsback to life of F_1) * and its Jacobian are calculated here in the rotated_back coordiantes*/ delta_P_K_seri = delta_stress_eval_F_1( mat_1_stress_eval, mat_2_stress_eval, F_coord, Eigen::Map(F_1_rot.data()), rotator, ratio); delta_P_seri = std::get<0>(delta_P_K_seri); delta_K_seri = std::get<1>(delta_P_K_seri); delta_P_eq = get_equation_strain_from_serial(delta_P_seri); delta_K_eq = get_equation_stiffness_from_serial(delta_K_seri); // std::cout << "residual _before: " << std::endl // << delta_P_eq << std::endl // << std::endl; // // solution loop: do { // updating variables: F_1_eq = F_1_new_eq; // solving for ∇(δK) * [x_n+1 - x_n] = -δP F_1_new_eq = F_1_eq - delta_K_eq.ldlt().solve(delta_P_eq); // auto F_1_new_temp = F_1_eq - // delta_K_eq.colPivHouseholderQr().solve(delta_P_eq); F_1_new_eq = F_1_eq // - delta_K_eq.colPivHouseholderQr().solve(delta_P_eq); F_1_new_seri = get_serial_strain_from_equation(F_1_new_eq); // updating F_1 before stress and stiffness evaluation F_1_rot = get_total_strain(F_1_new_seri, F_para); // stress and stiffenss evaluation in the rotated coordinates delta_P_K_seri = delta_stress_eval_F_1( mat_1_stress_eval, mat_2_stress_eval, F_coord, Eigen::Map(F_1_rot.data()), rotator, ratio); delta_P_seri = std::get<0>(delta_P_K_seri); delta_K_seri = std::get<1>(delta_P_K_seri); delta_P_eq = get_equation_strain_from_serial(delta_P_seri); delta_K_eq = get_equation_stiffness_from_serial(delta_K_seri); // computing variable used for loop termination criterion // energy norm of the residual - // std::cout << "iter : " << iter + 1 << std::endl << std::endl; - // std::cout << "residual : " << std::endl - // << delta_P_eq << std::endl - // << std::endl; - // std::cout << "delta_K_eq : " << std::endl - // << delta_K_eq << std::endl - // << std::endl; + std::cout << "iter : " << iter + 1 << std::endl << std::endl; + std::cout << "residual : " << std::endl + << delta_P_eq << std::endl + << std::endl; + std::cout << "delta_K_eq : " << std::endl + << delta_K_eq << std::endl + << std::endl; // std::cout << "F_new : " << std::endl // << F_1_new_eq - F_1_eq << std::endl << std::endl // << F_1_new_temp - F_1_eq << std::endl // << std::endl; del_energy = std::abs(delta_P_eq.dot(F_1_new_eq - F_1_eq)); + std::cout << "del_energy : " << del_energy << std::endl << std::endl; iter++; } while (del_energy > tol and iter < max_iter); // check if the loop has converged or not: if (iter == max_iter) { throw std::runtime_error("Error: the solver has not converged!!!!"); } // computing stress and stiffness of each layer according to it's // correspondent strin calculated in the loop F_1_rot = get_total_strain(F_1_new_seri, F_para); // storing the resultant strain in each layer in F_1 and F_2 auto F_2_new_seri = linear_eqs_seri(ratio, F_0_rot, F_1_new_seri); auto F_2_rot = get_total_strain(F_2_new_seri, F_para); F_1_mat = rotator.rotate(Eigen::Map(F_1_rot.data())); F_2_mat = rotator.rotate(Eigen::Map(F_2_rot.data())); P_K_1 = mat_1_stress_eval(F_1_mat); P_K_2 = mat_2_stress_eval(F_2_mat); P_1 = std::get<0>(P_K_1); P_1 = rotator.rotate_back(P_1.eval()); P_2 = std::get<0>(P_K_2); P_2 = rotator.rotate_back(P_2.eval()); K_1 = std::get<1>(P_K_1); K_1 = rotator.rotate_back(K_1.eval()); K_2 = std::get<1>(P_K_2); K_2 = rotator.rotate_back(K_2.eval()); // combine the computed strains and tangents to have the resultant // stress and tangent of the pixel ret_P = lam_P_combine(P_1, P_2, ratio); ret_K = lam_K_combine(K_1, K_2, ratio); return std::make_tuple(iter, del_energy, rotator.rotate(ret_P), rotator.rotate(ret_K)); } /* ---------------------------------------------------------------------- */ template struct LamHomogen; template struct LamHomogen; } // namespace muSpectre diff --git a/tests/split_test_laminate_solver.cc b/tests/split_test_laminate_solver.cc index 07dec6b..5cd1854 100644 --- a/tests/split_test_laminate_solver.cc +++ b/tests/split_test_laminate_solver.cc @@ -1,394 +1,394 @@ /** * @file test_laminate_solver.cc * * @author Till Junge * * @date 19 Oct 2018 * * @brief Tests for the large-strain, laminate homogenisation algorithm * * Copyright © 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 "materials/material_linear_elastic1.hh" #include "materials/material_linear_elastic2.hh" #include "materials/material_orthotropic.hh" #include "materials/laminate_homogenisation.hh" #include "tests.hh" #include "tests/test_goodies.hh" #include "common/field_collection.hh" #include "common/iterators.hh" #include "materials/materials_toolbox.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(laminate_homogenisation); template struct MaterialFixture { // constructor MaterialFixture(); }; // Material orthotropic fixture template struct MaterialFixture, Dim, c> { MaterialFixture() : mat("Name_aniso", aniso_inp_maker()) {} std::vector aniso_inp_maker() { std::vector aniso_inp{}; switch (Dim) { case (twoD): { aniso_inp = {c * 1.1, c * 2.1, c * 3.1, c * 4.1, c * 5.1, c * 6.1}; break; } case (threeD): { aniso_inp = {c * 1.1, c * 2.1, c * 3.1, c * 4.1, c * 5.1, c * 6.1, c * 7.1, c * 8.1, c * 9.1, c * 10.1, c * 11.1, c * 12.1, c * 13.1, c * 14.1, c * 15.1, c * 16.1, c * 17.1, c * 18.1, c * 19.1, c * 20.1, c * 21.1}; break; } default: throw std::runtime_error("Unknown dimension"); } return aniso_inp; } MaterialAnisotropic mat; }; // Material orthotropic fixture template struct MaterialFixture, Dim, c> { MaterialFixture() : mat("Name_orthro", ortho_inp_maker()) {} std::vector ortho_inp_maker() { std::vector ortho_inp{}; switch (Dim) { case (twoD): { ortho_inp = {c * 1.1, c * 2.1, c * 3.1, c * 4.1}; break; } case (threeD): { ortho_inp = {c * 1.1, c * 2.1, c * 3.1, c * 4.1, c * 5.1, c * 6.1, c * 7.1, c * 8.1, c * 9.1}; break; } default: throw std::runtime_error("Unknown dimension"); } return ortho_inp; } MaterialOrthotropic mat; }; // Material linear elastic fixture template struct MaterialFixture, Dim, c> { MaterialFixture() : mat("Name_LinElastic1", young_maker(), poisson_maker()) {} Real young_maker() { Real lambda{c * 2}, mu{c * 1.5}; Real young{mu * (3 * lambda + 2 * mu) / (lambda + mu)}; return young; } Real poisson_maker() { Real lambda{c * 2}, mu{c * 1.5}; Real poisson{lambda / (2 * (lambda + mu))}; return poisson; } MaterialLinearElastic1 mat; }; // Material pair fixture template struct MaterialPairFixture { using Vec_t = Eigen::Matrix; using Stiffness_t = T4Mat; using Serial_stiffness_t = Eigen::Matrix; using Parallel_stiffness_t = T4Mat; using Strain_t = Eigen::Matrix; using Serial_strain_t = Eigen::Matrix; using Parallel_strain_t = Eigen::Matrix; using Stress_t = Strain_t; using Serial_stress_t = Serial_strain_t; using Paralle_stress_t = Parallel_strain_t; using StrainMat_t = Eigen::Matrix; using StressMat_t = StrainMat_t; using T4MatMat_t = Eigen::Map>; using Output_t = std::tuple; using Function_t = std::function &)>; // constructor : MaterialPairFixture() { auto mat_fix_1 = std::unique_ptr>(); auto mat_fix_2 = std::unique_ptr>(); this->normal_vec = this->normal_vec / this->normal_vec.norm(); } constexpr Dim_t get_dim() { return Dim; } MaterialFixture mat_fix_1{}; MaterialFixture mat_fix_2{}; - Vec_t normal_vec{Vec_t::Random()}; - // Vec_t normal_vec{ Vec_t::UnitX()}; + // Vec_t normal_vec{Vec_t::Random()}; + Vec_t normal_vec{ (Vec_t::UnitY() + -3*Vec_t::UnitX()) /sqrt(10)}; // Real ratio{static_cast(std::rand() / (RAND_MAX))}; Real ratio = 0.5; static constexpr Dim_t fix_dim{Dim}; }; /** using fix_list = boost::mpl::list< MaterialPairFixture, MaterialLinearElastic1, threeD, 1, 1>, MaterialPairFixture, MaterialOrthotropic, threeD, 1, 1>, MaterialPairFixture, MaterialAnisotropic, threeD, 1, 1>, MaterialPairFixture, MaterialLinearElastic1, twoD, 1, 1>, MaterialPairFixture, MaterialOrthotropic, twoD, 1, 1>, MaterialPairFixture, MaterialAnisotropic, twoD, 1, 1>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(identical_material, Fix, fix_list, Fix) { auto & mat1{Fix::mat_fix_1.mat}; auto & mat2{Fix::mat_fix_2.mat}; using Strain_t = typename Fix::Strain_t; using StrainMat_t = typename Fix::StrainMat_t; StrainMat_t F; F = StrainMat_t::Random(); F = F.eval() + F.transpose().eval(); using Fucntion_t = typename Fix::Function_t; Fucntion_t mat1_evaluate_stress_func = [&mat1](const Eigen::Ref & strain) { return mat1.evaluate_stress_tangent(std::move(strain)); }; Fucntion_t mat2_evaluate_stress_func = [&mat2](const Eigen::Ref & strain) { return mat2.evaluate_stress_tangent(std::move(strain)); }; auto P_K_lam = LamHomogen::laminate_solver( Eigen::Map(F.data()), mat1_evaluate_stress_func, mat2_evaluate_stress_func, Fix::ratio, Fix::normal_vec, 1e-8, 1e5); auto P_lam = std::get<2>(P_K_lam); auto K_lam = std::get<3>(P_K_lam); auto P_K_ref = mat1_evaluate_stress_func(F); auto P_ref = std::get<0>(P_K_ref); auto K_ref = std::get<1>(P_K_ref); auto err_P{(P_lam - P_ref).norm()}; auto err_K{(K_lam - K_ref).norm()}; auto iters = std::get<0>(P_K_lam); auto del_energy = std::get<1>(P_K_lam); BOOST_CHECK_LT(err_P, tol); BOOST_CHECK_LT(err_K, tol); BOOST_CHECK_EQUAL(1, iters); BOOST_CHECK_LT(std::abs(del_energy), tol); }; **/ using fix_list3 = boost::mpl::list< /** MaterialPairFixture, MaterialLinearElastic1, threeD, 7, 3>, MaterialPairFixture, MaterialLinearElastic1, threeD, 3, 7>, MaterialPairFixture, MaterialLinearElastic1, twoD, 3, 1>, MaterialPairFixture, MaterialLinearElastic1, twoD, 1, 3>, MaterialPairFixture, MaterialOrthotropic, threeD, 7, 3>, MaterialPairFixture, MaterialOrthotropic, threeD, 3, 7>, MaterialPairFixture, MaterialOrthotropic, twoD, 3, 1>,**/ MaterialPairFixture, MaterialOrthotropic, twoD, 1, 2>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(different_material_PK_2, Fix, fix_list3, Fix) { auto & mat1{Fix::mat_fix_1.mat}; auto & mat2{Fix::mat_fix_2.mat}; using Strain_t = typename Fix::Strain_t; using StrainMat_t = typename Fix::StrainMat_t; StrainMat_t F; F = StrainMat_t::Zero(); F(0, 0) = 0.01; F(1, 0) = 0.005; - F(0, 1) = 0.005; + F(0, 1) = 0.007; F(1, 1) = 0.01; using Function_t = typename Fix::Function_t; Function_t mat1_evaluate_stress_func = [&mat1](const Eigen::Ref & strain) { return mat1.evaluate_stress_tangent(std::move(strain)); }; Function_t mat2_evaluate_stress_func = [&mat2](const Eigen::Ref & strain) { return mat2.evaluate_stress_tangent(std::move(strain)); }; - + std::cout << Fix::normal_vec << std::endl; auto solver_response = LamHomogen::laminate_solver( Eigen::Map(F.data()), mat1_evaluate_stress_func, mat2_evaluate_stress_func, Fix::ratio, Fix::normal_vec, 1e-9, 100); // auto P_lam = std::get<2> (P_K_lam); // auto K_lam = std::get<3>(P_K_lam); auto iters = std::get<0>(solver_response); auto del_energy = std::get<1>(solver_response); // auto P_K_1 = mat1_evaluate_stress_func(F); // auto P_1 = std::get<2> (P_K_1); // auto K_1 = std::get<3> (P_K_1); // auto P_K_2 = mat2_evaluate_stress_func(F); // auto P_2 = std::get<2> (P_K_2); // auto K_2 = std::get<3> (P_K_2); // auto K_ref = K_lam; // auto P_ref = Fix::get_mean_P(P_1, P_2, Fix::ratio); // auto err_K{(K_lam - K_ref).norm()}; // BOOST_CHECK_LT(err_P, tol); // BOOST_CHECK_LT(err_K, tol); - // BOOST_CHECK_EQUAL(1, iters); + BOOST_CHECK_EQUAL(1, iters); BOOST_CHECK_LT(0.0, del_energy); } /** using fix_list2 = boost::mpl::list , MaterialLinearElastic1, threeD, 7, 3>, MaterialPairFixture, MaterialLinearElastic1, threeD, 3, 7>, MaterialPairFixture, MaterialLinearElastic1, twoD, 3, 1>, MaterialPairFixture, MaterialLinearElastic1, twoD, 1, 3>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(different_material_PK_1, Fix, fix_list2, Fix) { auto & mat1{Fix::mat_fix_1.mat}; auto & mat2{Fix::mat_fix_2.mat}; using Strain_t = typename Fix::Strain_t; using StrainMat_t = typename Fix::StrainMat_t; StrainMat_t F; F = StrainMat_t::Zero(); F(0, 0) = 0.01; F(1, 0) = 0.005; F(0, 1) = 0.005; F(1, 1) = 0.01; using Function_t = typename Fix::Function_t; Dim_t Dim = Fix::get_dim(); Function_t mat1_evaluate_stress_func = [Dim, &mat1](const Eigen::Ref & strain){ // here the material calculates the stress in its own traits conversion // from strain to stress: auto mat_stress_tgt = mat1.evaluate_stress_tangent(std::move(strain)); // here we extract the type of strain and stress that mat1 deals with using traits = typename std::remove_reference_t::traits; const auto strain_measure{traits::strain_measure}; const auto stress_measure{traits::stress_measure}; auto ret_fun = [Dim, mat_stress_tgt, stress_measure, strain_measure] (const Eigen::Ref & strain){ auto && ret_P_K = MatTB::PK1_stress (strain, std::get<2>(mat_stress_tgt), std::get<3>(mat_stress_tgt)); return ret_P_K; }; return ret_fun(std::move(strain)); }; Function_t mat2_evaluate_stress_func = [Dim, &mat2](const Eigen::Ref & strain){ // here we extract the type of strain and stress that mat2 deals with using traits = typename std::remove_reference_t::traits; const auto strain_measure{traits::strain_measure}; const auto stress_measure{traits::stress_measure}; // here the material calculates the stress in its own traits conversion // from strain to stress: auto mat_stress_tgt = mat2.evaluate_stress_tangent(std::move(strain)); auto ret_fun = [Dim, mat_stress_tgt, stress_measure, strain_measure] (const Eigen::Ref & strain){ auto && ret_P_K = MatTB::PK1_stress (strain, std::get<2>(mat_stress_tgt), std::get<3>(mat_stress_tgt)); return ret_P_K; }; return ret_fun(std::move(strain)); }; auto P_K_lam = LamHomogen::laminate_solver(Eigen::Map(F.data()), mat1_evaluate_stress_func, mat2_evaluate_stress_func, Fix::ratio, Fix::normal_vec, 1e-9, 100); //auto P_lam = std::get<2> (P_K_lam); auto K_lam = std::get<3> (P_K_lam); //auto P_K_1 = mat1_evaluate_stress_func(F); //auto P_1 = std::get<2> (P_K_1); //auto K_1 = std::get<3> (P_K_1); //auto P_K_2 = mat2_evaluate_stress_func(F); //auto P_2 = std::get<2> (P_K_2); //auto K_2 = std::get<3> (P_K_2); auto K_ref = K_lam; //auto P_ref = Fix::get_mean_P(P_1, P_2, Fix::ratio); auto err_K{(K_lam - K_ref).norm()}; //BOOST_CHECK_LT(err_P, tol); BOOST_CHECK_LT(err_K, tol); } **/ BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_geometry.cc b/tests/test_geometry.cc index dd7987c..6d933ff 100644 --- a/tests/test_geometry.cc +++ b/tests/test_geometry.cc @@ -1,319 +1,324 @@ /** * file test_geometry.cc * * @author Till Junge * * @date 19 Apr 2018 * * @brief Tests for tensor rotations * * @section LICENSE * * Copyright © 2018 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 µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/geometry.hh" #include "tests.hh" #include "test_goodies.hh" #include "common/T4_map_proxy.hh" #include #include #include namespace muSpectre { enum class IsAligned { yes, no }; BOOST_AUTO_TEST_SUITE(geometry); /* ---------------------------------------------------------------------- */ template struct RotationFixture { static constexpr Dim_t Dim{Dim_}; using Vec_t = Eigen::Matrix; using Mat_t = Eigen::Matrix; using Ten_t = T4Mat; static constexpr Dim_t get_Dim() { return Dim_; } using Rot_t = RotatorBase; explicit RotationFixture(Rot_t rot_mat_inp) : rot_mat{rot_mat_inp}, rotator(rot_mat) {} Vec_t v{Vec_t::Random()}; Mat_t m{Mat_t::Random()}; Ten_t t{Ten_t::Random()}; Mat_t rot_mat; Rot_t rotator; testGoodies::RandRange rr{}; }; template struct RotationAngleFixture { static constexpr Dim_t Dim{Dim_}; using Parent = RotationFixture; using Vec_t = Eigen::Matrix; using Mat_t = Eigen::Matrix; using Ten_t = T4Mat; using Angles_t = Eigen::Matrix; using RotAng_t = RotatorAngle; static constexpr RotationOrder EulerOrder{Rot}; static constexpr Dim_t get_Dim() { return Dim_; } RotationAngleFixture() : rotator{euler} {} Vec_t v{Vec_t::Random()}; Mat_t m{Mat_t::Random()}; Ten_t t{Ten_t::Random()}; Angles_t euler{2 * pi * Angles_t::Random()}; RotatorAngle rotator; }; template struct RotationTwoVecFixture { static constexpr Dim_t Dim{Dim_}; using Parent = RotationFixture; using Vec_t = Eigen::Matrix; using Mat_t = Eigen::Matrix; using Ten_t = T4Mat; static constexpr Dim_t get_Dim() { return Dim_; } RotationTwoVecFixture() : vec_ref{this->ref_vec_maker()}, vec_des{this->des_vec_maker()}, rotator(vec_ref, vec_des) {} Vec_t ref_vec_maker() { Vec_t ret_vec{Vec_t::Random()}; return ret_vec / ret_vec.norm(); } Vec_t des_vec_maker() { if (is_aligned == IsAligned::yes) { return -this->vec_ref; } else { Vec_t ret_vec{Vec_t::Random()}; return ret_vec / ret_vec.norm(); } } Vec_t v{Vec_t::Random()}; Mat_t m{Mat_t::Random()}; Ten_t t{Ten_t::Random()}; Vec_t vec_ref{Vec_t::Random()}; Vec_t vec_des{Vec_t::Random()}; RotatorTwoVec rotator; }; + /* ---------------------------------------------------------------------- */ template struct RotationNormalFixture { static constexpr Dim_t Dim{Dim_}; using Parent = RotationFixture; using Vec_t = Eigen::Matrix; using Mat_t = Eigen::Matrix; using Ten_t = T4Mat; static constexpr Dim_t get_Dim() { return Dim_; } - RotationNormalFixture() : vec_norm{this->vec_maker()}, rotator(vec_norm) {} + + RotationNormalFixture() + : vec_norm{this->vec_maker()}, rotator(vec_norm) + {} Vec_t vec_maker() { if (is_aligned == IsAligned::yes) { return -Vec_t::UnitX(); } else { Vec_t ret_vec{Vec_t::Random()}; return ret_vec / ret_vec.norm(); } } + Vec_t v{Vec_t::Random()}; Mat_t m{Mat_t::Random()}; Ten_t t{Ten_t::Random()}; Vec_t vec_norm; RotatorNormal rotator; }; /* ---------------------------------------------------------------------- */ using fix_list = boost::mpl::list< RotationAngleFixture, RotationAngleFixture, RotationNormalFixture, RotationNormalFixture, RotationTwoVecFixture, RotationTwoVecFixture>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(rotation_test, Fix, fix_list, Fix) { using Vec_t = typename Fix::Vec_t; using Mat_t = typename Fix::Mat_t; using Ten_t = typename Fix::Ten_t; constexpr const Dim_t Dim{Fix::get_Dim()}; Vec_t & v{Fix::v}; Mat_t & m{Fix::m}; Ten_t & t{Fix::t}; const Mat_t & R{Fix::rotator.get_rot_mat()}; Vec_t v_ref{R * v}; Mat_t m_ref{R * m * R.transpose()}; Ten_t t_ref{Ten_t::Zero()}; for (int i = 0; i < Dim; ++i) { for (int a = 0; a < Dim; ++a) { for (int l = 0; l < Dim; ++l) { for (int b = 0; b < Dim; ++b) { for (int m = 0; m < Dim; ++m) { for (int n = 0; n < Dim; ++n) { for (int o = 0; o < Dim; ++o) { for (int p = 0; p < Dim; ++p) { get(t_ref, a, b, o, p) += R(a, i) * R(b, l) * get(t, i, l, m, n) * R(o, m) * R(p, n); } } } } } } } } Vec_t v_rotator(Fix::rotator.rotate(v)); Mat_t m_rotator(Fix::rotator.rotate(m)); Ten_t t_rotator(Fix::rotator.rotate(t)); auto v_error{(v_rotator - v_ref).norm() / v_ref.norm()}; BOOST_CHECK_LT(v_error, tol); auto m_error{(m_rotator - m_ref).norm() / m_ref.norm()}; BOOST_CHECK_LT(m_error, tol); auto t_error{(t_rotator - t_ref).norm() / t_ref.norm()}; BOOST_CHECK_LT(t_error, tol); if (t_error >= tol) { std::cout << "t4_reference:" << std::endl << t_ref << std::endl; std::cout << "t4_rotator:" << std::endl << t_rotator << std::endl; } Vec_t v_back{Fix::rotator.rotate_back(v_rotator)}; Mat_t m_back{Fix::rotator.rotate_back(m_rotator)}; Ten_t t_back{Fix::rotator.rotate_back(t_rotator)}; v_error = (v_back - v).norm() / v.norm(); BOOST_CHECK_LT(v_error, tol); m_error = (m_back - m).norm() / m.norm(); BOOST_CHECK_LT(m_error, tol); t_error = (t_back - t).norm() / t.norm(); BOOST_CHECK_LT(t_error, tol); } /* ---------------------------------------------------------------------- */ using threeD_list = boost::mpl::list< RotationAngleFixture>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(rotation_matrix_test, Fix, threeD_list, Fix) { using Mat_t = typename Fix::Mat_t; auto c{Eigen::cos(Fix::euler.array())}; Real c_1{c[0]}, c_2{c[1]}, c_3{c[2]}; auto s{Eigen::sin(Fix::euler.array())}; Real s_1{s[0]}, s_2{s[1]}, s_3{s[2]}; Mat_t rot_ref; switch (Fix::EulerOrder) { case RotationOrder::ZXYTaitBryan: { rot_ref << c_1 * c_3 - s_1 * s_2 * s_3, -c_2 * s_1, c_1 * s_3 + c_3 * s_1 * s_2, c_3 * s_1 + c_1 * s_2 * s_3, c_1 * c_2, s_1 * s_3 - c_1 * c_3 * s_2, -c_2 * s_3, s_2, c_2 * c_3; break; } default: { BOOST_CHECK(false); break; } } auto err{(rot_ref - Fix::rotator.get_rot_mat()).norm()}; BOOST_CHECK_LT(err, tol); if (not(err < tol)) { std::cout << "Reference:" << std::endl << rot_ref << std::endl; std::cout << "Rotator:" << std::endl << Fix::rotator.get_rot_mat() << std::endl; } } /* ---------------------------------------------------------------------- */ using twovec_list = boost::mpl::list, RotationTwoVecFixture, RotationTwoVecFixture, RotationTwoVecFixture>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(rotation_twovec_test, Fix, twovec_list, Fix) { using Vec_t = typename Fix::Vec_t; Vec_t vec_ref{Fix::vec_ref}; Vec_t vec_des{Fix::vec_des}; Vec_t vec_res{Fix::rotator.rotate(vec_ref)}; Vec_t vec_back{Fix::rotator.rotate_back(vec_res)}; auto err_f{(vec_res - vec_des).norm()}; BOOST_CHECK_LT(err_f, tol); if (err_f >= tol) { std::cout << "Destination:" << std::endl << vec_des << std::endl; std::cout << "Rotated:" << std::endl << vec_res << std::endl; } auto err_b{(vec_back - vec_ref).norm()}; BOOST_CHECK_LT(err_b, tol); if (err_b >= tol) { std::cout << "Refrence:" << std::endl << vec_ref << std::endl; std::cout << "Rotated Back:" << std::endl << vec_back << std::endl; } } /* ---------------------------------------------------------------------- */ using normal_list = boost::mpl::list, RotationNormalFixture, RotationNormalFixture, RotationNormalFixture>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(rotation_normal_test, Fix, normal_list, Fix) { using Vec_t = typename Fix::Vec_t; Vec_t vec_ref{Fix::vec_norm}; Vec_t vec_des{Vec_t::UnitX()}; Vec_t vec_res{Fix::rotator.rotate_back(vec_ref)}; Vec_t vec_back{Fix::rotator.rotate(vec_res)}; auto err_f{(vec_res - vec_des).norm()}; BOOST_CHECK_LT(err_f, tol); if (err_f >= tol) { std::cout << "Destination:" << std::endl << vec_des << std::endl; std::cout << "Rotated:" << std::endl << vec_res << std::endl; } auto err_b{(vec_back - vec_ref).norm()}; BOOST_CHECK_LT(err_b, tol); if (err_b >= tol) { std::cout << "Refrence:" << std::endl << vec_ref << std::endl; std::cout << "Rotated Back:" << std::endl << vec_back << std::endl; } } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_materials_toolbox.cc b/tests/test_materials_toolbox.cc index 97885a8..92fb176 100644 --- a/tests/test_materials_toolbox.cc +++ b/tests/test_materials_toolbox.cc @@ -1,371 +1,375 @@ /** * @file test_materials_toolbox.cc * * @author Till Junge * * @date 05 Nov 2017 * * @brief Tests for the materials toolbox * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include "tests.hh" #include "materials/materials_toolbox.hh" #include "materials/stress_transformations_default_case.hh" #include "materials/stress_transformations_PK1_impl.hh" #include "materials/stress_transformations_PK2_impl.hh" #include "materials/stress_transformations.hh" #include "common/T4_map_proxy.hh" #include "common/tensor_algebra.hh" #include "tests/test_goodies.hh" #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(materials_toolbox) BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_strain_conversion, Fix, testGoodies::dimlist, Fix) { constexpr Dim_t dim{Fix::dim}; using T2 = Eigen::Matrix; T2 F{(T2::Random() - .5 * T2::Ones()) / 20 + T2::Identity()}; // checking Green-Lagrange T2 Eref = .5 * (F.transpose() * F - T2::Identity()); T2 E_tb = MatTB::convert_strain( Eigen::Map>(F.data())); Real error = (Eref - E_tb).norm(); BOOST_CHECK_LT(error, tol); // checking Left Cauchy-Green Eref = F * F.transpose(); E_tb = MatTB::convert_strain(F); error = (Eref - E_tb).norm(); BOOST_CHECK_LT(error, tol); // checking Right Cauchy-Green Eref = F.transpose() * F; E_tb = MatTB::convert_strain(F); error = (Eref - E_tb).norm(); BOOST_CHECK_LT(error, tol); // checking Hencky (logarithmic) strain Eref = F.transpose() * F; Eigen::SelfAdjointEigenSolver EigSolv(Eref); Eref.setZero(); for (size_t i{0}; i < dim; ++i) { auto && vec = EigSolv.eigenvectors().col(i); auto && val = EigSolv.eigenvalues()(i); Eref += .5 * std::log(val) * vec * vec.transpose(); } E_tb = MatTB::convert_strain(F); error = (Eref - E_tb).norm(); BOOST_CHECK_LT(error, tol); auto F_tb = MatTB::convert_strain( F); error = (F - F_tb).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(dumb_tensor_mult_test, Fix, testGoodies::dimlist, Fix) { constexpr Dim_t dim{Fix::dim}; using T4 = T4Mat; T4 A, B, R1, R2; A.setRandom(); B.setRandom(); R1 = A * B; R2.setZero(); for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t a = 0; a < dim; ++a) { for (Dim_t b = 0; b < dim; ++b) { for (Dim_t k = 0; k < dim; ++k) { for (Dim_t l = 0; l < dim; ++l) { get(R2, i, j, k, l) += get(A, i, j, a, b) * get(B, a, b, k, l); } } } } } } auto error{(R1 - R2).norm()}; BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_PK1_stress, Fix, testGoodies::dimlist, Fix) { using Matrices::Tens2_t; using Matrices::Tens4_t; using Matrices::tensmult; constexpr Dim_t dim{Fix::dim}; using T2 = Eigen::Matrix; using T4 = T4Mat; testGoodies::RandRange rng; T2 F = T2::Identity() * 2; // F.setRandom(); T2 E_tb = MatTB::convert_strain( Eigen::Map>(F.data())); Real lambda = 3; // rng.randval(1, 2); Real mu = 4; // rng.randval(1,2); T4 J = Matrices::Itrac(); T2 I = Matrices::I2(); T4 I4 = Matrices::Isymm(); T4 C = lambda * J + 2 * mu * I4; T2 S = Matrices::tensmult(C, E_tb); T2 Sref = lambda * E_tb.trace() * I + 2 * mu * E_tb; auto error{(Sref - S).norm()}; BOOST_CHECK_LT(error, tol); T4 K = Matrices::outer_under(I, S) + Matrices::outer_under(F, I) * C * Matrices::outer_under(F.transpose(), I); // See Curnier, 2000, "Méthodes numériques en mécanique des solides", p 252 T4 Kref; Real Fkrkr = (F.array() * F.array()).sum(); T2 Fkmkn = F.transpose() * F; T2 Fisjs = F * F.transpose(); Kref.setZero(); for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t m = 0; m < dim; ++m) { for (Dim_t n = 0; n < dim; ++n) { get(Kref, i, m, j, n) = (lambda * ((Fkrkr - dim) / 2 * I(i, j) * I(m, n) + F(i, m) * F(j, n)) + mu * (I(i, j) * Fkmkn(m, n) + Fisjs(i, j) * I(m, n) - I(i, j) * I(m, n) + F(i, n) * F(j, m))); } } } } error = (Kref - K).norm(); BOOST_CHECK_LT(error, tol); T2 P = MatTB::PK1_stress( F, S); T2 Pref = F * S; error = (P - Pref).norm(); BOOST_CHECK_LT(error, tol); auto && stress_tgt = MatTB::PK1_stress( F, S, C); T2 P_t = std::move(std::get<0>(stress_tgt)); T4 K_t = std::move(std::get<1>(stress_tgt)); error = (P_t - Pref).norm(); BOOST_CHECK_LT(error, tol); error = (K_t - Kref).norm(); BOOST_CHECK_LT(error, tol); auto && stress_tgt_trivial = MatTB::PK1_stress(F, P, K); T2 P_u = std::move(std::get<0>(stress_tgt_trivial)); T4 K_u = std::move(std::get<1>(stress_tgt_trivial)); error = (P_u - Pref).norm(); BOOST_CHECK_LT(error, tol); error = (K_u - Kref).norm(); BOOST_CHECK_LT(error, tol); T2 P_g; T4 K_g; std::tie(P_g, K_g) = testGoodies::objective_hooke_explicit(lambda, mu, F); error = (P_g - Pref).norm(); BOOST_CHECK_LT(error, tol); error = (K_g - Kref).norm(); BOOST_CHECK_LT(error, tol); } BOOST_AUTO_TEST_CASE(elastic_modulus_conversions) { // define original input constexpr Real E{123.456}; constexpr Real nu{.3}; using MatTB::convert_elastic_modulus; // derived values constexpr Real K{ convert_elastic_modulus(E, nu)}; constexpr Real lambda{ convert_elastic_modulus(E, nu)}; constexpr Real mu{ convert_elastic_modulus(E, nu)}; // recover original inputs Real comp = convert_elastic_modulus(K, mu); Real err = E - comp; BOOST_CHECK_LT(err, tol); comp = convert_elastic_modulus(K, mu); err = nu - comp; BOOST_CHECK_LT(err, tol); comp = convert_elastic_modulus(lambda, mu); err = E - comp; BOOST_CHECK_LT(err, tol); // check inversion resistance Real compA = convert_elastic_modulus(K, mu); Real compB = convert_elastic_modulus(mu, K); BOOST_CHECK_EQUAL(compA, compB); // check trivial self-returning comp = convert_elastic_modulus(K, mu); BOOST_CHECK_EQUAL(K, comp); comp = convert_elastic_modulus(K, mu); BOOST_CHECK_EQUAL(mu, comp); // check alternative calculation of computed values comp = convert_elastic_modulus( K, mu); // alternative for "Shear" BOOST_CHECK_LE(std::abs((comp - lambda) / lambda), tol); } template struct FiniteDifferencesHolder { constexpr static FiniteDiff value{FinDiff}; }; using FinDiffList = boost::mpl::list, FiniteDifferencesHolder, FiniteDifferencesHolder>; - /* ---------------------------------------------------------------------- */ - BOOST_FIXTURE_TEST_CASE_TEMPLATE(numerical_tangent_test, Fix, FinDiffList, - Fix) { - constexpr Dim_t Dim{twoD}; - using T4_t = T4Mat; - using T2_t = Eigen::Matrix; - - bool verbose{false}; - - T4_t Q{}; - Q << 1., 2., 0., 0., 0., 1.66666667, 0., 0., 0., 0., 2.33333333, 0., 0., 0., - 0., 3.; - if (verbose) { - std::cout << Q << std::endl << std::endl; - } - - T2_t B{}; - B << 2., 3.33333333, 2.66666667, 4.; - if (verbose) { - std::cout << B << std::endl << std::endl; - } - - auto fun = [&](const T2_t & x) -> T2_t { - using cmap_t = Eigen::Map>; - using map_t = Eigen::Map>; - cmap_t x_vec{x.data()}; - T2_t ret_val{}; - map_t(ret_val.data()) = Q * x_vec + cmap_t(B.data()); - return ret_val; - }; - - T2_t temp_res = fun(T2_t::Ones()); - if (verbose) { - std::cout << temp_res << std::endl << std::endl; - } - - T4_t numerical_tangent{MatTB::compute_numerical_tangent( - fun, T2_t::Ones(), 1e-2)}; - - if (verbose) { - std::cout << numerical_tangent << std::endl << std::endl; - } - - Real error{(numerical_tangent - Q).norm() / Q.norm()}; - - BOOST_CHECK_LT(error, tol); - if (not(error < tol)) { - switch (Fix::value) { - case FiniteDiff::backward: { - std::cout << "backward difference: " << std::endl; - break; - } - case FiniteDiff::forward: { - std::cout << "forward difference: " << std::endl; - break; - } - case FiniteDiff::centred: { - std::cout << "centered difference: " << std::endl; - break; - } - } - - std::cout << "error = " << error << std::endl; - std::cout << "numerical tangent:\n" << numerical_tangent << std::endl; - std::cout << "reference:\n" << Q << std::endl; - } - } + /* --------------------------------------------------- + ------------------- */ + // // BOOST_FIXTURE_TEST_CASE_TEMPLATE(numerical_tangent_test, Fix, + // FinDiffList, + // // Fix) { + // // constexpr Dim_t Dim{twoD}; + // // using T4_t = T4Mat; + // // using T2_t = Eigen::Matrix; + + // // bool verbose{false}; + + // // T4_t Q{}; + // // Q << 1., 2., 0., 0., 0., 1.66666667, 0., 0., 0., 0., + // // 2.33333333, 0., 0., 0.,0., 3.; + // // if (verbose) { + // // std::cout << Q << std::endl << std::endl; + // // } + + // // T2_t B{}; + // // B << 2., 3.33333333, 2.66666667, 4.; + // // if (verbose) { + // // std::cout << B << std::endl << std::endl; + // // } + + // // auto fun = [&](const T2_t & x) -> T2_t { + // // using cmap_t = Eigen::Map>; + // // using map_t = Eigen::Map>; + // // cmap_t x_vec{x.data()}; + // // T2_t ret_val{}; + // // map_t(ret_val.data()) = Q * x_vec + cmap_t(B.data()); + // // return ret_val; + // // }; + + // // T2_t temp_res = fun(T2_t::Ones()); + // // if (verbose) { + // // std::cout << temp_res << std::endl << std::endl; + // // } + + // // T4_t numerical_tangent{MatTB::compute_numerical_tangent( + // // fun, T2_t::Ones(), 1e-2)}; + + // // if (verbose) { + // // std::cout << numerical_tangent << std::endl << std::endl; + // // } + + // // Real error{(numerical_tangent - Q).norm() / Q.norm()}; + + // // BOOST_CHECK_LT(error, tol); + // // if (not(error < tol)) { + // // switch (Fix::value) { + // // case FiniteDiff::backward: { + // // std::cout << "backward difference: " << std::endl; + // // break; + // // } + // // case FiniteDiff::forward: { + // // std::cout << "forward difference: " << std::endl; + // // break; + // // } + // // case FiniteDiff::centred: { + // // std::cout << "centered difference: " << std::endl; + // // break; + // // } + // // } + + // // std::cout << "error = " << error << std::endl; + // // std::cout << "numerical tangent:\n" << numerical_tangent << + // std::endl; + // // std::cout << "reference:\n" << Q << std::endl; + // // } + // // } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre