diff --git a/src/materials/laminate_homogenisation.cc b/src/materials/laminate_homogenisation.cc index e927c16..2ea68a5 100644 --- a/src/materials/laminate_homogenisation.cc +++ b/src/materials/laminate_homogenisation.cc @@ -1,583 +1,598 @@ /** * @file lamminate_hemogenization.hh * * @author Ali Falsafi * * @date 28 Sep 2018 * * @brief * * 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 "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_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; } /* ---------------------------------------------------------------------- */ // It does not to be especialised because all functions used init are // especialised. this two functions compute the strai in the second // layer of laminate material from strain in the first layer and total strain template auto LamHomogen::linear_eqs (Real ratio, const Eigen::Ref& F_lam, const Eigen::Ref& F_1) ->Strain_t{ Strain_t ret_F; Serial_strain_t F_lam_seri = get_serial_strain(F_lam); Serial_strain_t F_1_seri = get_serial_strain(F_1); Serial_strain_t F_2_seri{}; F_2_seri = (F_lam_seri - ratio * F_1_seri) / (1-ratio); Parallel_strain_t F_lam_para = get_parallel_strain(F_lam); Strain_t F_2 = get_total_strain(F_2_seri, F_lam_para); return F_2; } template auto LamHomogen:: linear_eqs_seri (Real ratio, const Eigen::Ref & F_lam, const Eigen::Ref & F_1) ->Serial_strain_t{ Serial_strain_t F_2_seri; Serial_strain_t F_lam_seri = get_serial_strain(F_lam); F_2_seri = (F_lam_seri - ratio * F_1) / (1-ratio); return F_2_seri; } /* ---------------------------------------------------------------------- */ 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) ->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(del_P_mat); auto del_P_seri = get_serial_strain(Eigen::Map (del_P_mat.data()) ); Stiffness_t del_K; del_K = std::get<1>(P_K_1) - std::get<1>(P_K_2); del_K = rotator.rotate(del_K); auto del_K_seri = get_serial_stiffness(del_K); return std::make_tuple(del_P_seri, del_K_seri); } /* ---------------------------------------------------------------------- */ 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, const Real ratio) ->StressMat_t{ StressMat_t P_tmp{}; //return P_tmp ; 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()); } /* ---------------------------------------------------------------------- */ 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}; RotatorNormal rotator(normal_vec); Strain_t F_1{F_coord},F_2{F_coord}; Parallel_strain_t F_para{}; F_para = get_parallel_strain(F_coord); Serial_strain_t F_1_new_seri{get_serial_strain(F_coord)}; + Serial_strain_t F_2_new_seri{linear_eqs_seri(ratio, F_coord, F_1_new_seri)}; Serial_strain_t F_1_seri {F_1_new_seri}; - Serial_strain_t F_2_seri {F_1_new_seri}; + Serial_strain_t F_2_seri {F_2_new_seri}; std::tuple delta_P_K_seri; Serial_stress_t delta_P_seri{}; Serial_stiffness_t delta_K_seri{}; 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 and its Jacobian is calculated here */ delta_P_K_seri = delta_stress_eval(mat_1_stress_eval, mat_2_stress_eval, Eigen::Map (F_1.data()), Eigen::Map (F_2.data()), rotator); delta_P_seri = std::get<0> (delta_P_K_seri); delta_K_seri = std::get<1> (delta_P_K_seri); - /*Serial_stiffness_t A{Serial_stiffness_t::Identity()};// + 1e-1*Serial_stiffness_t::Random()}; - Serial_strain_t b{Serial_strain_t::Random()}; - auto f{[&A,&b](auto x){ return A*x + b;}}; - auto grad_f{[&A](auto x){ return A;}}; - delta_K_seri = grad_f(F_1_seri); - delta_P_seri = f(F_1_seri);*/ + /*Serial_stiffness_t A1{Serial_stiffness_t::Identity()};// + 1e-1*Serial_stiffness_t::Random()}; + auto A2 = 2 * A1; + Serial_strain_t b1{Serial_strain_t::Zero()};//{Serial_strain_t::Random()}; + auto b2 = 2 * b1; + auto f1{[&A1,&b1](auto x){ return A1*x + b1;}}; + auto grad_f1{[&A1](auto x){ return A1;}}; + + auto f2{[&A2,&b2](auto x){ return A2*x + b2;}}; + auto grad_f2{[&A2](auto x){ return A2;}};*/ + + //delta_K_seri = grad_f2(F_1_seri)-grad_f1(F_1_seri); + //delta_P_seri = f2(F_1_seri) - f1(F_1_seri); //std::cout<<" b :" << b.transpose() << std::endl; std::cout<<" del_P_0 :" << delta_P_seri.transpose() << std::endl; std::cout<<" F_1_seri_0 :" << F_1_seri.transpose() << std::endl; + do { // updating variables: - F_1_seri = F_1_new_seri; - F_2_seri = linear_eqs_seri(ratio, F_coord, F_1_seri); - F_1 = get_total_strain(F_1_seri, F_para); - F_2 = get_total_strain(F_2_seri, F_para); + F_1_seri = F_1_new_seri; + F_2_seri = F_2_new_seri; + std::cout<<" del_P :" << delta_P_seri.transpose() << std::endl; + std::cout<<" del_K_old :" << std::endl<< delta_K_seri.transpose() << std::endl; //solving for ∇(delP) * [x_n+1 - x_n] = -delP - F_1_new_seri = F_1_seri - delta_K_seri.ldlt().solve(delta_P_seri); + F_1_new_seri = F_1_seri - delta_K_seri.ldlt().solve(delta_P_seri); + 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); delta_P_K_seri = delta_stress_eval(mat_1_stress_eval, mat_2_stress_eval, Eigen::Map (F_1.data()), Eigen::Map (F_2.data()), rotator); delta_P_seri = std::get<0> (delta_P_K_seri); delta_K_seri = std::get<1> (delta_P_K_seri); - /* delta_K_seri = grad_f(F_1_new_seri); - delta_P_seri = f(F_1_new_seri);*/ + /*delta_K_seri = grad_f2(F_1_seri) - grad_f1(F_1_seri); + delta_P_seri = f2(F_1_seri) - f1(F_1_seri);*/ + // computing variable used for loop termination criterion // energy norm of the residual del_energy = std::abs(delta_P_seri.dot(F_1_new_seri - F_1_seri)); - std::cout<<" step :" << iter << std::endl; + std::cout<< std::endl<<" step :" << iter << std::endl; + std::cout<<" F_1 :" << F_1.transpose() << std::endl; + std::cout<<" F_2 :" << F_2.transpose() << std::endl; + std::cout<<" F_old :" << F_1_seri.transpose() << std::endl; + std::cout<<" F_new :" << F_1_new_seri.transpose() << std::endl; std::cout<<" residual :" << delta_P_seri.transpose() << std::endl; - std::cout<<" del_K :" << delta_K_seri.transpose() << std::endl; + std::cout<<" del_K :" << std::endl<< delta_K_seri.transpose() << std::endl; std::cout<< "energy_norm :" << del_energy << std::endl; iter ++; } while(del_energy > tol && iter < max_iter); - std::cout << "ieterations number : "<< iter << std::endl; - std::cout<<" del_P :" << delta_P_seri.transpose() << std::endl; - std::cout<<" del_K :" << delta_K_seri.transpose() << std::endl; - std::cout<<" F_1_seri_1 :" << F_1_seri.transpose() << std::endl; + std::cout << "iterations number : "<< iter << std::endl; + std::cout<<" final residual :" << delta_P_seri.transpose() << std::endl; + std::cout<<" final_F_seri_1 :" << F_1_seri.transpose() << std::endl; // storing the reuktant strain in each layer in F_1 and F_2 F_1 = get_total_strain(F_1_seri, F_para); F_2 = get_total_strain(F_2_seri, F_para); // computing stress and stiffness of each layer according to it's // correspondent strin calculated in the loop P_K_1 = mat_1_stress_eval(Eigen::Map (F_1.data())); P_K_2 = mat_2_stress_eval(Eigen::Map (F_2.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); // combine the computed strains and tangents to have the resultant // stress and tangent of the pixel P_1 = rotator.rotate(P_1); P_2 = rotator.rotate(P_2); ret_P = lam_P_combine(P_1, P_2, ratio); K_1 = rotator.rotate(K_1); K_2 = rotator.rotate(K_2); ret_K = lam_K_combine(K_1, K_2, ratio); return std::make_tuple(rotator.rotate_back(ret_P), rotator.rotate_back(ret_K)); //(ret_P,ret_K); } /* ---------------------------------------------------------------------- */ template struct LamHomogen; template struct LamHomogen; } //muSpectre diff --git a/tests/split_test_laminate_solver.cc b/tests/split_test_laminate_solver.cc index 7290b3c..9a735f5 100644 --- a/tests/split_test_laminate_solver.cc +++ b/tests/split_test_laminate_solver.cc @@ -1,385 +1,385 @@ /** * @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(); }; 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; }; 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; }; 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; }; 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 StrainMap_t = Eigen::Matrix; using StressMap_t = StrainMap_t; using T4MatMap_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>(); } constexpr Dim_t get_dim () {return Dim;} Stiffness_t get_mean_K(Stiffness_t K_1, Stiffness_t K_2, Real ratio){ Stiffness_t ret_K{Stiffness_t::Zero()}; for (int i = 0 ; i < Dim; ++i) { for (int j = 0; j < Dim; ++j) { for (int k = 0; k < Dim; ++k) { for (int l = 0; l < Dim; ++l) { auto& retk = get(ret_K, i, j, k, l); const auto& k1 = get(K_1, i, j, k, l); const auto& k2 = get(K_2, i, j, k, l); if (i == 0 or j == 0 or k == 0 or l == 0) { retk = 1 / ( ratio / k1 + (1-ratio) / k2 ) ; } else { retk = ratio * k1 + (1-ratio) * k2 ; } } } } } return ret_K; } StressMap_t get_mean_P(StressMap_t P_1, StressMap_t P_2, Real ratio){ StressMap_t ret_P{StressMap_t::Zero()}; for (int i = 0 ; i < Dim; ++i) { for (int j = 0; j < Dim; ++j) { if (i == 0 or j == 0){ ret_P(i,j) = 1 / (1/(ratio * P_1(i, j)) + 1/((1-ratio) * P_2(i, j)) ) ; } else { ret_P(i,j) = ratio * P_1(i,j) + ((1-ratio)*P_2(i, j)); } } } std::cout << "P_1 = "<< std::endl << P_1 << std::endl; std::cout << "P_2 = " << std::endl<< P_2 << std::endl; return ret_P; } MaterialFixture mat_fix_1{}; MaterialFixture mat_fix_2{}; //Vec_t normal_vec{Vec_t::Random()}; Vec_t normal_vec{Vec_t::UnitX()}; //Real ratio {(double) std::rand() / (RAND_MAX)}; Real ratio = 0.5; static constexpr Dim_t fix_dim{Dim}; }; /* 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>>; 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 StrainMap_t = typename Fix::StrainMap_t; StrainMap_t F; F = StrainMap_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)); }; 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<0> (P_K_lam); auto K_lam = std::get<1> (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()}; BOOST_CHECK_LT(err_P, tol); BOOST_CHECK_LT(err_K, tol); }; */ using fix_list2 = boost::mpl::list , MaterialLinearElastic1, threeD, 7, 3>, MaterialPairFixture, MaterialLinearElastic1, threeD, 3, 7>>;/*, MaterialPairFixture, MaterialLinearElastic1, twoD, 1, 3>>;*/ BOOST_FIXTURE_TEST_CASE_TEMPLATE(different_material, 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 StrainMap_t = typename Fix::StrainMap_t; StrainMap_t F; //F = StrainMap_t::Random(); F<< 0.01, 0, 0, 0 , 0, 0, - 0 , 0, 0.01; + 0 , 0, 0.0; Dim_t Dim = Fix::get_dim(); using Function_t = typename Fix::Function_t; // using Output_t = typename Fix::Output_t; 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)); 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)); 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-6, + 1e-9, 100); //auto P_lam = std::get<0> (P_K_lam); auto K_lam = 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_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 K_ref = Fix::get_mean_K(K_1, K_2, Fix::ratio); //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(); } //muSpectre