Page MenuHomec4science
Files F63274447

split_test_laminate_solver.cc
No OneTemporary
Actions

File Metadata

Created
Sat, May 18, 23:36

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();
};
// Material orthotropic fixture
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;
};
// Material orthotropic fixture
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;
};
// Material linear elastic fixture
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;
};
// Material pair fixture
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();
}
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{static_cast<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, 1, 1>,
MaterialPairFixture<MaterialAnisotropic<threeD, threeD>,
MaterialAnisotropic<threeD, threeD>, threeD, 1, 1>,
MaterialPairFixture<MaterialLinearElastic1<twoD, twoD>,
MaterialLinearElastic1<twoD, twoD>, twoD, 1, 1>,
MaterialPairFixture<MaterialOrthotropic<twoD, twoD>,
MaterialOrthotropic<twoD, twoD>, twoD, 1, 1>,
MaterialPairFixture<MaterialAnisotropic<twoD, twoD>,
MaterialAnisotropic<twoD, twoD>, 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<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<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_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<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<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<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<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<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);
}
**/
using fix_list3 = boost::mpl::
/*MaterialPairFixture<MaterialLinearElastic1<threeD, threeD>,
MaterialLinearElastic1<threeD, threeD>,
threeD, 7, 3>,
MaterialPairFixture<MaterialLinearElastic1<threeD, threeD>,
MaterialLinearElastic1<threeD, threeD>,
threeD, 3, 7>>;*/
list<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<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<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_LT(0.0, del_energy);
}
BOOST_AUTO_TEST_SUITE_END();
} // namespace muSpectre

Event Timeline

c4science · Help