Homec4science
Diffusion µSpectre (master) Browse

muSpectre/tests/split_test_laminate_solver.cc
1af5123b9fbamaster

Last Change
Raw File

split_test_laminate_solver.cc
View Options

/**
* @file test_laminate_solver.cc
*
* @author Till Junge <ali.falsafi@epfl.ch>
*
* @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 <type_traits>
#include <boost/mpl/list.hpp>
#include <boost/range/combine.hpp>
#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 <class Mat_t, Dim_t Dim, Dim_t c>
struct MaterialFixture{
//constructor
MaterialFixture();
};
template<Dim_t Dim, Dim_t c>
struct MaterialFixture<MaterialAnisotropic<Dim, Dim>, Dim, c>{
MaterialFixture():
mat("Name_aniso", aniso_inp_maker()){}
std::vector<Real> aniso_inp_maker(){
std::vector<Real> 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<Dim, Dim> mat;
};
template<Dim_t Dim, Dim_t c>
struct MaterialFixture<MaterialOrthotropic<Dim, Dim>, Dim, c>{
MaterialFixture():
mat("Name_orthro", ortho_inp_maker()){}
std::vector<Real> ortho_inp_maker(){
std::vector<Real> 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<Dim, Dim> mat;
};
template<Dim_t Dim, Dim_t c>
struct MaterialFixture<MaterialLinearElastic1<Dim, Dim>, 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<Dim, Dim> mat;
};
template <class Mat1_t, class Mat2_t, Dim_t Dim, Dim_t c1, Dim_t c2>
struct MaterialPairFixture
{
using Vec_t = Eigen::Matrix<Real, Dim, 1>;
using Stiffness_t = T4Mat<Real, Dim>;
using Serial_stiffness_t = Eigen::Matrix<Real, 2*Dim-1, 2*Dim-1>;
using Parallel_stiffness_t = T4Mat<Real, (Dim-1)>;
using Strain_t = Eigen::Matrix<Real, Dim*Dim, 1>;
using Serial_strain_t = Eigen::Matrix<Real, 2*Dim-1, 1>;
using Parallel_strain_t = Eigen::Matrix<Real, (Dim-1)*(Dim-1), 1>;
using Stress_t = Strain_t;
using Serial_stress_t = Serial_strain_t;
using Paralle_stress_t = Parallel_strain_t;
using StrainMap_t = Eigen::Matrix<Real, Dim, Dim>;
using StressMap_t = StrainMap_t;
using T4MatMap_t = Eigen::Map<T4Mat<Real, Dim>>;
using Output_t = std::tuple<StressMap_t,
Stiffness_t>;
using Function_t = std::function<Output_t(const Eigen::Ref<StrainMap_t> &)>;
// constructor :
MaterialPairFixture(){
auto mat_fix_1 = std::unique_ptr<MaterialFixture<Mat1_t, Dim, c1>>();
auto mat_fix_2 = std::unique_ptr<MaterialFixture<Mat2_t, Dim, c2>>();
}
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<Mat1_t, Dim, c1> mat_fix_1{};
MaterialFixture<Mat2_t, Dim, c2> 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
<MaterialPairFixture<MaterialLinearElastic1<threeD, threeD>,
MaterialLinearElastic1<threeD, threeD>,
threeD, 1, 1>,
MaterialPairFixture<MaterialOrthotropic<threeD, threeD>,
MaterialOrthotropic<threeD, threeD>,
threeD, 2, 2>,
MaterialPairFixture<MaterialAnisotropic<threeD, threeD>,
MaterialAnisotropic<threeD, threeD>,
threeD, 2, 2>,
MaterialPairFixture<MaterialLinearElastic1<twoD, twoD>,
MaterialLinearElastic1<twoD, twoD>,
twoD, 2, 2>,
MaterialPairFixture<MaterialOrthotropic<twoD, twoD>,
MaterialOrthotropic<twoD, twoD>,
twoD, 1, 1>,
MaterialPairFixture<MaterialAnisotropic<twoD, twoD>,
MaterialAnisotropic<twoD, twoD>,
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<StrainMap_t> & strain){
return mat1.evaluate_stress_tangent(std::move(strain));
};
Fucntion_t mat2_evaluate_stress_func =
[&mat2](const Eigen::Ref<StrainMap_t> & strain){
return mat2.evaluate_stress_tangent(std::move(strain));
};
auto P_K_lam =
LamHomogen<Fix::fix_dim>::laminate_solver(Eigen::Map<Strain_t>(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
<MaterialPairFixture<MaterialLinearElastic1<threeD, threeD>,
MaterialLinearElastic1<threeD, threeD>,
threeD, 7, 3>,
MaterialPairFixture<MaterialLinearElastic1<threeD, threeD>,
MaterialLinearElastic1<threeD, threeD>,
threeD, 3, 7>>;/*,
MaterialPairFixture<MaterialLinearElastic1<twoD, twoD>,
MaterialLinearElastic1<twoD, twoD>,
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;
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<StrainMap_t> & 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<decltype(mat1)>::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<StrainMap_t> & strain){
auto && ret_P_K = MatTB::PK1_stress<stress_measure, strain_measure>
(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<StrainMap_t> & strain){
// here we extract the type of strain and stress that mat2 deals with
using traits = typename std::remove_reference_t<decltype(mat2)>::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<StrainMap_t> & strain){
auto && ret_P_K = MatTB::PK1_stress<stress_measure, strain_measure>
(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<Fix::fix_dim>::laminate_solver(Eigen::Map<Strain_t>(F.data()),
mat1_evaluate_stress_func,
mat2_evaluate_stress_func,
Fix::ratio,
Fix::normal_vec,
1e-6,
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
c4science · Help