Page MenuHomec4science

split_test_laminate_solver.cc
No OneTemporary

File Metadata

Created
Mon, Nov 4, 16:44

split_test_laminate_solver.cc

/**
* @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 StrainMat_t = Eigen::Matrix<Real, Dim, Dim>;
using StressMat_t = StrainMat_t;
using T4MatMat_t = Eigen::Map<T4Mat<Real, Dim>>;
using Output_t = std::tuple<StressMat_t,
Stiffness_t>;
using Function_t = std::function<Output_t(const Eigen::Ref<StrainMat_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>>();
this->normal_vec = this->normal_vec / this->normal_vec.norm();
//std::cout << this-> normal_vec;
}
constexpr Dim_t get_dim () {return Dim;}
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 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<StrainMat_t> & strain){
return mat1.evaluate_stress_tangent(std::move(strain));
};
Fucntion_t mat2_evaluate_stress_func =
[&mat2](const Eigen::Ref<StrainMat_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, 3, 1>,
MaterialPairFixture<MaterialLinearElastic1<twoD, twoD>,
MaterialLinearElastic1<twoD, twoD>,
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<StrainMat_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<StrainMat_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<StrainMat_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<StrainMat_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-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 = 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::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, 3, 1>,
MaterialPairFixture<MaterialLinearElastic1<twoD, twoD>,
MaterialLinearElastic1<twoD, twoD>,
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<StrainMat_t>& strain){
return mat1.evaluate_stress_tangent(std::move(strain));};
Function_t mat2_evaluate_stress_func = [&mat2]
(const Eigen::Ref<StrainMat_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-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 = 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();
} //muSpectre

Event Timeline