diff --git a/src/materials/laminate_homogenisation.cc b/src/materials/laminate_homogenisation.cc index 6bed6c1..1b0dbdd 100644 --- a/src/materials/laminate_homogenisation.cc +++ b/src/materials/laminate_homogenisation.cc @@ -1,628 +1,596 @@ /** * @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 (unsigned int i = 0; i < Dim; ++i) { for (unsigned int j = 0; j < Dim; ++j) { ret_K(i, j) = Ks(i, j); } for (unsigned int j = 1; j < Dim; ++j) { ret_K(i, j) += Ks(i, j + (Dim - 1)); } } return ret_K; } /* ---------------------------------------------------------------------- */ // 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) -> 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) / (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; del_K = (std::get<1>(P_K_1) - (-ratio / (1 - ratio)) * std::get<1>(P_K_2)); del_K = rotator.rotate_back(del_K).eval(); 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(Eigen::Ref P_1, - Eigen::Ref P_2, + auto LamHomogen::lam_P_combine(const Eigen::Ref& P_1, + const Eigen::Ref& P_2, const Real ratio) -> StressMat_t { - StressMat_t P_tmp{}; - // return P_tmp ; - Serial_stress_t P_1_seri = + /*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_P = get_total_strain(P_1_seri, P_para); - return Eigen::Map(ret_P.data()); + 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 { + -> 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); - // ///////////temporary : - // StrainMat_t aaaa = - StrainMat_t::Identity(); - // rotator.set_rot_mat(aaaa); - // /////////// - Strain_t F_1{F_coord}, F_2{F_coord}; StrainMat_t F_mat{rotator.rotate_back(Eigen::Map(F_1.data()))}; 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); - // monitoring - /*std::cout << std::endl << "F_1 Initial : " << F_1.transpose() << - * std::endl;*/ - 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 rotate_backd coordiantes - */ + /* (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); - // monitoring : - /* std::cout << "residual before the loop : " << delta_P_eq.transpose() - << 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); 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 del_energy = std::abs(delta_P_eq.dot(F_1_new_eq - F_1_eq)); - // monitoring : - /* std::cout << " step : " << iter + 1 << std::endl - << " del_energy " - << " of step " << iter + 1 << " is :" << del_energy << std::endl - << " residual " - << " of step " << iter + 1 << " is :" << delta_P_eq.transpose() - << std::endl - << " del_K " - << " of step " << iter + 1 << " is :" << std::endl - << delta_K_eq << std::endl;*/ iter++; } while (del_energy > tol && iter < max_iter); // check if the loop has lead in convergence 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); - auto F_2_rot = linear_eqs_seri(ratio, F_0_rot, F_1_new_seri); - F_mat = rotator.rotate(Eigen::Map(F_1_rot.data())); - F_1 = Eigen::Map(F_mat.data()); - F_mat = rotator.rotate(Eigen::Map(F_2_rot.data())); - F_2 = Eigen::Map(F_mat.data()); - P_K_1 = mat_1_stress_eval(Eigen::Map(F_1.data())); - P_K_2 = mat_2_stress_eval(Eigen::Map(F_2.data())); + + // 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); + P_K_1 = mat_1_stress_eval(Eigen::Map(F_1_rot.data())); + P_K_2 = mat_2_stress_eval(Eigen::Map(F_2_rot.data())); P_1 = std::get<0>(P_K_1); P_2 = std::get<0>(P_K_2); K_1 = std::get<1>(P_K_1); K_2 = std::get<1>(P_K_2); - // storing the resultant strain in each layer in F_1 and F_2 - auto F_2_new_seri = linear_eqs_seri(ratio, F_coord, F_1_new_seri); - F_1 = get_total_strain(F_1_new_seri, F_para); - F_2 = get_total_strain(F_2_new_seri, F_para); - // monitoring : - /* std::cout << "steps : " << iter << std::endl; - std::cout << "F_1 Final : " << F_1.transpose() << std::endl; - std::cout << "F_2 Final : " << F_2.transpose() << std::endl; */ - // 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(rotator.rotate(ret_P), rotator.rotate(ret_K)); + 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/src/materials/laminate_homogenisation.hh b/src/materials/laminate_homogenisation.hh index 3a5859a..836a1f8 100644 --- a/src/materials/laminate_homogenisation.hh +++ b/src/materials/laminate_homogenisation.hh @@ -1,202 +1,202 @@ /** * @file lamminate_hemogenization.hh * * @author Ali Falsafi * * @date 28 Sep 2018 * * @brief Laminatehomogenisation enables one to obtain the resulting stress * and stiffness tensors of a laminate pixel that is consisted of two * materialswith a certain normal vector of their interface plane. * note that it is supposed to be used in static way. so it does note * any data member. It is merely a collection of functions used to * calculate effective stress and stiffness. * * * 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. */ #ifndef SRC_MATERIALS_LAMINATE_HOMOGENISATION_HH_ #define SRC_MATERIALS_LAMINATE_HOMOGENISATION_HH_ #include "common/geometry.hh" #include "common/common.hh" #include "common/field_map.hh" namespace muSpectre { // TODO(faslafi): to template with the formulation as well template struct LamHomogen { //! typedefs for data handled by this interface using Vec_t = Eigen::Matrix; using Stiffness_t = T4Mat; using Serial_stiffness_t = Eigen::Matrix; using Parallel_stiffness_t = T4Mat; using Equation_stiffness_t = Eigen::Matrix; using Strain_t = Eigen::Matrix; using Serial_strain_t = Eigen::Matrix; using Parallel_strain_t = Eigen::Matrix; using Equation_strain_t = Eigen::Matrix; using Stress_t = Strain_t; using Serial_stress_t = Serial_strain_t; using Parallel_stress_t = Parallel_strain_t; using Equation_stress_t = Equation_strain_t; using StrainMat_t = Eigen::Matrix; using StressMat_t = StrainMat_t; using Function_t = std::function( const Eigen::Ref &)>; // these two functions return the indexes in column major strain and stress // tensors that the behavior is either serial or parallel concerning // combining the laminate layers static constexpr auto get_serial_indexes() -> std::array; static constexpr auto get_parallel_indexes() -> std::array; static constexpr auto get_if_serial() -> std::array; // these functions return the parts of a stress or strain tensor that // behave either serial or parallel in a laminate structure // they are used to obtain a certian part of tensor which is needed by the // solver used in this struct to calculate the effective stress and // stiffness or compose the complete tensor from its sub sets. // obtain the serial part of stress or strain tensor static auto get_serial_strain(const Eigen::Ref & F) -> Serial_strain_t; // obtain the serial part of stress or strain tensor static auto get_equation_strain_from_serial(const Eigen::Ref & F) -> Equation_strain_t; static auto get_serial_strain_from_equation(const Eigen::Ref & Fs) -> Serial_strain_t; // obtain the parallel part of stress or strain tensor static auto get_parallel_strain(const Eigen::Ref & F) -> Parallel_strain_t; // compose the complete strain or stress tensor from // its serial and parallel parts static auto get_total_strain(const Eigen::Ref & F_seri, const Eigen::Ref & F_para) -> Strain_t; // obtain the serial part of stifness tensor static auto get_serial_stiffness(const Eigen::Ref & K) -> Serial_stiffness_t; // obtain the parallel part of stifness tensor static auto get_parallel_stiffness(const Eigen::Ref & K) -> Parallel_stiffness_t; static auto get_equation_stiffness_from_serial( const Eigen::Ref & Ks) -> Equation_stiffness_t; // compose the complete stiffness tensor from // its serial and parallel parts static auto get_total_stiffness(const Eigen::Ref & K_seri, const Eigen::Ref & K_para) -> Stiffness_t; /** * this functions calculate the strain in the second materials considering * the fact that: * 1. In the directions that the laminate layers are in series the weighted * average of the strains should be equal to the total strain of the * laminate * 2. In the directions that the laminate layers are in parallel the strain * is constant and equal in all of the layers. */ static auto linear_eqs(Real ratio, const Eigen::Ref & F_0, const Eigen::Ref & F_1) -> Strain_t; static auto linear_eqs_seri(Real ratio, const Eigen::Ref & F_0, const Eigen::Ref & F_1) -> Serial_strain_t; /** * the objective in hemgenisation of a single laminate pixel is equating the * stress in the serial directions so the difference of stress between their * layers should tend to zero. this function return the stress difference * and the difference of Stiffness matrices which is used as the Jacobian in * the solution process */ static auto 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; static auto 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; /** * the following function s claculates the energy dimension error of the * solution. it will beused in each step of the solution to detremine the * relevant difference that implementation of that step has had on * convergence to the solution. */ inline static auto del_energy_eval(const Real del_F_norm, const Real delta_P_norm) -> Real; /** * This functions calculate the resultant stress and tangent matrices * according to the computed F_1 and F_2 from the solver. * */ - static auto lam_P_combine(Eigen::Ref P_1, - Eigen::Ref P_2, const Real ratio) - -> StressMat_t; + static auto lam_P_combine(const Eigen::Ref & P_1, + const Eigen::Ref & P_2, + const Real ratio) -> StressMat_t; static auto lam_K_combine(const Eigen::Ref & K_1, const Eigen::Ref & K_2, const Real ratio) -> Stiffness_t; /** * This is the main solver function that might be called staically from * an external file. this will return the resultant stress and stiffness * tensor according to interanl "equilibrium" of the lamiante. * The inputs are : * 1- global Strain * 2- stress calculation function of the layer 1 * 3- stress calculation function of the layer 2 * 4- the ratio of the first material in the laminate sturucture of the * pixel 5- the normal vector of the interface of two layers 6- the * tolerance error for the internal solution of the laminate pixel 7- the * maximum iterations for the internal solution of the laminate pixel */ static auto 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 = 1e-6, Dim_t max_iter = 1000) - -> std::tuple; + -> std::tuple; }; // LamHomogen } // namespace muSpectre #endif // SRC_MATERIALS_LAMINATE_HOMOGENISATION_HH_ diff --git a/tests/split_test_laminate_solver.cc b/tests/split_test_laminate_solver.cc index 73a0259..cc076e2 100644 --- a/tests/split_test_laminate_solver.cc +++ b/tests/split_test_laminate_solver.cc @@ -1,386 +1,395 @@ /** * @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()}; - Real ratio{static_cast(std::rand() / (RAND_MAX))}; - // Real ratio = 0.5; + // 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, twoD, 2, 2>>; /* - using fix_list = boost::mpl::list - , - MaterialLinearElastic1, - threeD, 1, 1>, - MaterialPairFixture, - MaterialOrthotropic, - threeD, 2, 2>, - MaterialPairFixture, - MaterialAnisotropic, - threeD, 2, 2>, - MaterialPairFixture, - MaterialLinearElastic1, - twoD, 2, 2>, - MaterialPairFixture, - MaterialOrthotropic, - twoD, 1, 1>, - MaterialPairFixture, - MaterialAnisotropic, - twoD, 2, 2>>; + , + MaterialLinearElastic1, + threeD, 1, 1>, + MaterialPairFixture, + MaterialOrthotropic, + threeD, 2, 2>, + MaterialPairFixture, + MaterialAnisotropic, + threeD, 2, 2>, + MaterialPairFixture, + MaterialLinearElastic1, + twoD, 2, 2>, + MaterialPairFixture, + MaterialOrthotropic, + twoD, 1, 1>, + MaterialPairFixture, + MaterialAnisotropic, + twoD, 2, 2>>;*/ 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(); 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)); - }; + [&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)); - }; + [&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<0> (P_K_lam); - auto K_lam = std::get<1> (P_K_lam); + 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 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_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<0>(mat_stress_tgt), std::get<1>(mat_stress_tgt)); + (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<0>(mat_stress_tgt), std::get<1>(mat_stress_tgt)); + (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<0> (P_K_lam); - auto K_lam = std::get<1> (P_K_lam); + //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<0> (P_K_1); - //auto K_1 = std::get<1> (P_K_1); + //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<0> (P_K_2); - //auto K_2 = std::get<1> (P_K_2); + //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); } **/ using fix_list3 = boost::mpl:: /*MaterialPairFixture, MaterialLinearElastic1, threeD, 7, 3>, MaterialPairFixture, MaterialLinearElastic1, threeD, 3, 7>>;*/ list, MaterialLinearElastic1, twoD, 3, 1>, MaterialPairFixture, MaterialLinearElastic1, twoD, 1, 3>>; 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(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)); }; 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<0> (P_K_lam); - auto K_lam = std::get<1>(P_K_lam); + // auto P_lam = std::get<2> (P_K_lam); + auto K_lam = std::get<3>(P_K_lam); + auto iters = std::get<0>(P_K_lam); + auto del_energy = std::get<1>(P_K_lam); // auto P_K_1 = mat1_evaluate_stress_func(F); - // auto P_1 = std::get<0> (P_K_1); - // auto K_1 = std::get<1> (P_K_1); + // 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<0> (P_K_2); - // auto K_2 = std::get<1> (P_K_2); + // 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_LT(0.0, del_energy); } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre