Homec4science
Diffusion µSpectre prototype implementation (add_eigenstrain) Browse

muSpectre_proto/tests/test_material_hyperelastic.cc
8353160ae1bdadd_eigenstrain

Last Change
Raw File

test_material_hyperelastic.cc
View Options

/**
* file test_material_hyperelastic.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 01 May 2017
*
* @brief Test the hyper-elastic material used in Section 4 of
* de Geus et al. / Comput. Methods in Appl. Mech. Engrg.
* 318, (2017), 412–430 https://doi.org/10.1016/j.cma.2016.12.032
*
* @section LICENCE
*
* Copyright (C) 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 "materials/material_hyper_elastic.hh"
#include <boost/test/unit_test.hpp>
#include <boost/mpl/list.hpp>
#include <stdexcept>
#include "tests.hh"
namespace muSpectre {
// params from "material_hyper_elastic.ipynb";
const Real E = 1.9058231542461002;
const Real nu = 0.2993067590987868;
template<Dim_t s_dim, Dim_t m_dim>
struct mat_fixture: public MaterialHyperElastic<s_dim, m_dim>{
mat_fixture():
MaterialHyperElastic<s_dim, m_dim>("material", E, nu){
};
const static Dim_t DimS = s_dim;
const static Dim_t DimM = m_dim;
/** inputs with expected outputs computed with
* "material_hyper_elastic.ipynb" in folder "helpers"*/
using parent = MaterialHyperElastic<s_dim, m_dim>;
using SecondTens = typename parent::SecondTens;
using FourthVoigt = Eigen::Matrix<Real, m_dim*m_dim, m_dim*m_dim>;
const static SecondTens F, P;
const static FourthVoigt K;
};
//----------------------------------------------------------------------------//
template <> const mat_fixture<2, 2>::SecondTens mat_fixture<2, 2>::F =
(mat_fixture<2, 2>::SecondTens() <<
1.0, 0.1,
0.2, 1.0).finished();
template <> const mat_fixture<2, 2>::SecondTens mat_fixture<2, 2>::P =
(mat_fixture<2, 2>::SecondTens()<<
0.07868216666666673, 0.22348781666666673,
0.2313560333333334, 0.07868216666666672).finished();
template <> const mat_fixture<2, 2>::FourthVoigt mat_fixture<2, 2>::K =
(mat_fixture<2, 2>::FourthVoigt()<<
2.624580833333334, 1.1084346666666667, 0.40273666666666674, 0.5854533333333334,
1.1084346666666667, 2.6245808333333334, 0.40273666666666674, 0.5854533333333334,
0.5854533333333334, 0.5854533333333334, 0.7552753333333334, 0.8925028333333335,
0.40273666666666674, 0.40273666666666674, 0.7936838333333334, 0.7552753333333334).finished();
//----------------------------------------------------------------------------//
template <> const mat_fixture<2, 3>::SecondTens mat_fixture<2, 3>::F =
(mat_fixture<2, 3>::SecondTens()<<
1.0, 0.1, 0.0,
0.2, 1.0, 0.0,
0.0, 0.0, 1.0).finished();
template <> const mat_fixture<2, 3>::SecondTens mat_fixture<2, 3>::P =
(mat_fixture<2, 3>::SecondTens()<<
0.07868216666666673, 0.22348781666666673, 0.0,
0.2313560333333334, 0.07868216666666672, 0.0,
0.0, 0.0, 0.027344166666666694).finished();
template <> const mat_fixture<2, 3>::FourthVoigt mat_fixture<2, 3>::K =
(mat_fixture<2, 3>::FourthVoigt()<<2.624580833333334, 1.1084346666666667, 1.0937666666666668, 0.0, 0.0, 0.40273666666666674, 0.0, 0.0, 0.5854533333333334,
1.1084346666666667, 2.6245808333333334, 1.0937666666666668, 0.0, 0.0, 0.40273666666666674, 0.0, 0.0, 0.5854533333333334,
1.0937666666666668, 1.0937666666666668, 2.5879108333333334, 0.0, 0.0, 0.10937666666666668, 0.0, 0.0, 0.21875333333333336,
0.0, 0.0, 0.0, 0.7334, 0.07334, 0.0, 0.7680781666666667, 0.22002000000000005, 0.0,
0.0, 0.0, 0.0, 0.14668, 0.7334, 0.0, 0.22002000000000005, 0.7900801666666668, 0.0,
0.5854533333333334, 0.5854533333333334, 0.21875333333333336, 0.0, 0.0, 0.7552753333333334, 0.0, 0.0, 0.8925028333333335,
0.0, 0.0, 0.0, 0.7900801666666668, 0.22002, 0.0, 0.7334, 0.14668, 0.0,
0.0, 0.0, 0.0, 0.22002, 0.7680781666666667, 0.0, 0.07334, 0.7334, 0.0,
0.40273666666666674, 0.40273666666666674, 0.10937666666666668, 0.0, 0.0, 0.7936838333333334, 0.0, 0.0, 0.7552753333333334).finished();
//----------------------------------------------------------------------------//
template <> const mat_fixture<3, 3>::SecondTens mat_fixture<3, 3>::F =
(mat_fixture<3, 3>::SecondTens()<<
1.0, 0.1, 0.3,
0.2, 1.0, 0.4,
0.5, 0.6, 1.0
).finished();
template <> const mat_fixture<3, 3>::SecondTens mat_fixture<3, 3>::P =
(mat_fixture<3, 3>::SecondTens()<<
0.9479714333333337, 0.7435627833333334, 0.92523635,
0.8402667666666668, 1.1591906333333337, 1.1568859333333334,
1.264590916666667, 1.4368351000000001, 1.4569510333333335).finished();
template <> const mat_fixture<3, 3>::FourthVoigt mat_fixture<3, 3>::K =
(mat_fixture<3,3>::FourthVoigt() <<
3.3442565, 1.1084346666666667, 1.2037766666666667, 0.4815106666666667, 1.193542, 0.6227566666666668, 0.69293, 1.5443073333333335, 0.6734613333333335,
1.1084346666666667, 3.4762685000000006, 1.2697826666666667, 1.4862686666666667, 0.35746600000000006, 0.4907446666666667, 1.90304, 0.6348913333333334, 0.8054733333333335,
1.2037766666666667, 1.2697826666666667, 3.6889545000000004, 1.5376066666666668, 1.178874, 0.2413886666666667, 1.851702, 1.5589753333333336, 0.3654333333333334,
0.69293, 1.90304, 1.851702, 0.9959040000000001, 0.270218, 0.7403540000000001, 2.6075758333333336, 0.9881900000000001, 0.49795200000000006,
1.5443073333333335, 0.6348913333333334, 1.5589753333333336, 0.3654333333333334, 0.8974650000000001, 0.4947283333333334, 0.9881900000000001, 2.3479155, 0.9894566666666668,
0.6734613333333335, 0.8054733333333335, 0.3654333333333334, 0.7915653333333335, 0.358986, 0.7552753333333334, 0.49795200000000006, 0.9894566666666668, 1.6635165000000003,
0.4815106666666667, 1.4862686666666667, 1.5376066666666668, 1.8534405000000003, 0.5272880000000001, 0.2637706666666667, 0.9959040000000001, 0.3654333333333334, 0.7915653333333335,
1.193542, 0.35746600000000006, 1.178874, 0.5272880000000001, 1.6521988333333337, 0.810217, 0.270218, 0.8974650000000001, 0.358986,
0.6227566666666668, 0.4907446666666667, 0.2413886666666667, 0.2637706666666667, 0.810217, 1.5940335000000003, 0.7403540000000001, 0.4947283333333334, 0.7552753333333334).finished();
typedef boost::mpl::list<mat_fixture<2, 2>,
mat_fixture<2, 3>,
mat_fixture<3, 3>> test_types;
BOOST_FIXTURE_TEST_CASE_TEMPLATE(material_hyper_elastic_test, F, test_types, F) {
BOOST_CHECK_CLOSE(F::lambda, E*nu/((1+nu)*(1-2*nu)), tol);
BOOST_CHECK_CLOSE(F::mu, E/2/(1+nu), tol);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_elastic_tensor, F, test_types, F) {
Eigen::Matrix<Real, vsize(F::DimM), vsize(F::DimM)> Cv;
F::VoigtConv::fourth_to_voigt(F::C, Cv);
auto & E = F::Young;
auto & nu = F::Poisson;
auto prefactor = E/(1+nu)/(1-2*nu);
const auto diff{vsize(F::DimM) - F::DimM};
for (Dim_t i = 0; i < F::DimM; ++i) {
BOOST_CHECK_CLOSE(Cv(i, i ), prefactor*(1-nu), 1e-12);
for (Dim_t j = i + 1; j < F::DimM; ++j ) {
BOOST_CHECK_CLOSE(Cv(i, j), prefactor*nu, 1e-12);
BOOST_CHECK_CLOSE(Cv(j, i), prefactor*nu, 1e-12);
}
}
for (Dim_t i = 0; i < diff; ++i) {
BOOST_CHECK_CLOSE(Cv(i+F::DimM, i+F::DimM), prefactor*(1-2*nu)/2, 1e-12);
for (Dim_t j = i + 1; j < diff; ++j ) {
BOOST_CHECK_CLOSE(Cv(i+F::DimM, j+F::DimM), 0, 1e-12);
}
for (Dim_t j = 0; j < diff; ++j) {
BOOST_CHECK_CLOSE(Cv(i+F::DimM, j ), 0, 1e-12);
BOOST_CHECK_CLOSE(Cv(i, j+F::DimM), 0, 1e-12);
}
}
}
//----------------------------------------------------------------------------//
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constitutive_law, F, test_types, F) {
//using SecondArray = typename MaterialHyperElastic<F::DimS, F::DimM>::SecondArray;
//using FourthArray = typename MaterialHyperElastic<F::DimS, F::DimM>::FourthArray;
using SecondTens = typename MaterialHyperElastic<F::DimS, F::DimM>::SecondTens;
using FourthTens = typename MaterialHyperElastic<F::DimS, F::DimM>::FourthTens;
//auto F_t = getMap(add_dim<decltype(F::F), F::DimS>(F::F));
//auto K_t = Eigen::TensorMap<FourthArray>(const_cast<Real*>(F::K.data()), F::DimM, F::DimM, F::DimM, F::DimM, 1, 1);
//std::cout << F_t;
//std::cout << K_t;
this->add_pixel(Ccoord_t<F::DimS>{0});
this->initialize();
std::array<Dim_t, F::DimS> sizes;
for (auto & s: sizes) {
s=1;
}
auto F_t = F::get_second_array_map(F::F.data(), sizes);
//std::cout << "hello: " << std::endl << F_t << std::endl;
SecondTens P_;
FourthTens K_;
auto P_t = F::get_second_array_map(P_.data(), sizes);
auto K_t = F::get_fourth_array_map(K_.data(), sizes);
for (Dim_t i = 0; i < F::DimM; ++i) {
for (Dim_t j = 0; j < F::DimM; ++j) {
P_(i, j) = 10*(i+1)+j+1;
}
}
F::compute_First_Piola_Kirchhoff_stress(F_t, P_t, K_t);
auto K_v = VoigtConversion<F::DimM>::template fourth_to_voigt<FourthTens, false>(K_.shuffle(std::array<int, 4>{1,0,2,3}));
Real P_error = (P_-F::P).norm();
if (P_error >= tol) {
std::cout << "First Piola-Kirhoff stress wrong:" << std::endl
<< "P ref = " << std::endl << F::P << std::endl
<< "P comp = " << std::endl << P_ << std::endl
<< "P diff = " << std::endl << F::P-P_ << std::endl<< std::endl;
}
BOOST_CHECK_LT(P_error, tol);
Real K_error = (K_v - F::K).norm();
if (K_error >= tol) {
std::cout << "Tangent wrong: "<< std::endl
<< "K ref = " << std::endl << F::K <<std::endl
<< "K comp = " << std::endl << K_v << std::endl
<< "K diff = " << std::endl << K_v - F::K << std::endl << std::endl;
}
BOOST_CHECK_LT(K_error, tol);
}
} // muSpectre
c4science · Help