Homec4science
Diffusion µSpectre (feat/repair_intersection_octree) Browse

muSpectre/tests/test_material_orthotropic.cc
5be13b02c658feat/repair_intersectio…

Last Change
Raw File

test_material_orthotropic.cc
View Options

/**
* @file test_material_orhtotropic.cc
*
* @author Ali Falsafi <ali.falsafi@epfl.ch>
*
* @date 20 Jul 2018
*
* @brief Tests for the standard Newton-Raphson + Conjugate Gradient solver
*
* 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 "tests.hh"
#include "solver/solvers.hh"
#include "solver/solver_cg.hh"
#include "solver/solver_eigen.hh"
#include "fft/fftw_engine.hh"
#include "fft/projection_finite_strain_fast.hh"
#include "materials/material_linear_elastic1.hh"
#include "materials/material_orthotropic.hh"
#include "common/iterators.hh"
#include "common/ccoord_operations.hh"
#include "common/common.hh"
#include "cell/cell_factory.hh"
#include "solver/deprecated_solvers.hh"
#include "solver/deprecated_solver_cg.hh"
#include "solver/deprecated_solver_cg_eigen.hh"
#include <boost/mpl/list.hpp>
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(anisotropic_material_tests);
BOOST_AUTO_TEST_CASE(orthotropic_twoD) {
constexpr Dim_t dim{twoD};
constexpr Ccoord_t<dim> resolutions{5, 5};
constexpr Rcoord_t<dim> lengths{5, 5};
auto fft_ptr_non{std::make_unique<FFTWEngine<dim>>(resolutions, ipow(dim, 2))};
auto proj_ptr_non{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr_non), lengths)};
CellBase<dim, dim> sys_non(std::move(proj_ptr_non));
auto fft_ptr_lin{std::make_unique<FFTWEngine<dim>>(resolutions, ipow(dim, 2))};
auto proj_ptr_lin{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr_lin), lengths)};
CellBase<dim, dim> sys_lin(std::move(proj_ptr_lin));
using Mat_t_non = MaterialOrthotropic<dim, dim>;
using Mat_t_lin = MaterialLinearElastic1<dim, dim>;
const Real Young{1e10}, Poisson{.3};
const int con {2} ;
const Real lambda{Young*Poisson/((1+Poisson)*(1-2*Poisson))};
const Real mu{Young/(2*(1+Poisson))};
auto& Material_soft_lin = Mat_t_lin::make(sys_lin, "soft_lin", Young, Poisson);
auto& Material_hard_lin = Mat_t_lin::make(sys_lin, "hard_lin", con * Young, Poisson);
std::vector<Real> input_soft;
input_soft.push_back(lambda+2*mu);
input_soft.push_back(lambda);
input_soft.push_back(lambda+2*mu);
input_soft.push_back(mu);
std::vector<Real> input_hard;
input_hard.push_back(con *(lambda+2*mu));
input_hard.push_back(con * lambda);
input_hard.push_back(con *(lambda+2*mu));
input_hard.push_back(con * mu);
auto& Material_soft = Mat_t_non::make(sys_non, "soft", input_soft);
auto& Material_hard = Mat_t_non::make(sys_non, "hard", input_hard);
for (auto && tup: akantu::enumerate(sys_non)) {
auto && pixel = std::get<1>(tup);
if (std::get<0>(tup) == 0) {
Material_soft.add_pixel(pixel);
Material_soft_lin.add_pixel(pixel);
} else {
Material_hard.add_pixel(pixel);
Material_hard_lin.add_pixel(pixel);
}
}
sys_lin.initialise();
sys_non.initialise();
Grad_t<dim> delF0;
delF0 << 1e-4, 5e-5, 5e-5, 0;//, 0, 0, 0, 0, 0;
constexpr Real cg_tol{1e-8}, newton_tol{1e-5};
constexpr Uint maxiter{CcoordOps::get_size(resolutions)*ipow(dim, secondOrder)*10};
constexpr bool verbose{false};
GradIncrements<dim> grads; grads.push_back(delF0);
DeprecatedSolverCG<dim> cg_lin{sys_lin, cg_tol, maxiter, bool(verbose)};
Eigen::ArrayXXd res2{deprecated_newton_cg(sys_lin, grads, cg_lin, newton_tol, verbose)[0].grad};
DeprecatedSolverCG<dim> cg_non{sys_non, cg_tol, maxiter, bool(verbose)};
Eigen::ArrayXXd res1{deprecated_newton_cg(sys_non, grads, cg_non, newton_tol, verbose)[0].grad};
BOOST_CHECK_LE(abs(res1-res2).mean(), cg_tol);
}
BOOST_AUTO_TEST_CASE(orthotropic_threeD) {
constexpr Dim_t dim{threeD};
constexpr Ccoord_t<dim> resolutions{5, 5, 5};
constexpr Rcoord_t<dim> lengths{5, 5, 5};
auto fft_ptr_non{std::make_unique<FFTWEngine<dim>>(resolutions, ipow(dim, 2))};
auto proj_ptr_non{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr_non), lengths)};
CellBase<dim, dim> sys_non(std::move(proj_ptr_non));
auto fft_ptr_lin{std::make_unique<FFTWEngine<dim>>(resolutions, ipow(dim, 2))};
auto proj_ptr_lin{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr_lin), lengths)};
CellBase<dim, dim> sys_lin(std::move(proj_ptr_lin));
using Mat_t_non = MaterialOrthotropic<dim, dim>;
using Mat_t_lin = MaterialLinearElastic1<dim, dim>;
const Real Young{1e10}, Poisson{.3};
//const Real E{1.0}, n{0.33};
const int con {2} ;
const Real lambda{Young*Poisson/((1+Poisson)*(1-2*Poisson))};
const Real mu{Young/(2*(1+Poisson))};
auto& Material_soft_lin = Mat_t_lin::make(sys_lin, "soft_lin", Young, Poisson);
auto& Material_hard_lin = Mat_t_lin::make(sys_lin, "hard_lin", con * Young, Poisson);
std::vector<Real> input_soft;
input_soft.push_back(lambda+2*mu);
input_soft.push_back(lambda);
input_soft.push_back(lambda);
input_soft.push_back(lambda+2*mu);
input_soft.push_back(lambda);
input_soft.push_back(lambda+2*mu);
input_soft.push_back(mu);
input_soft.push_back(mu);
input_soft.push_back(mu);
std::vector<Real> input_hard;
input_hard.push_back(con *(lambda+2*mu));
input_hard.push_back(con * lambda);
input_hard.push_back(con * lambda);
input_hard.push_back(con *(lambda+2*mu));
input_hard.push_back(con * lambda);
input_hard.push_back(con *(lambda+2*mu));
input_hard.push_back(con * mu);
input_hard.push_back(con * mu);
input_hard.push_back(con * mu);
auto& Material_soft = Mat_t_non::make(sys_non, "soft", input_soft);
auto& Material_hard = Mat_t_non::make(sys_non, "hard", input_hard);
for (auto && tup: akantu::enumerate(sys_non)) {
auto && pixel = std::get<1>(tup);
if (std::get<0>(tup) == 0) {
Material_soft.add_pixel(pixel);
Material_soft_lin.add_pixel(pixel);
} else {
Material_hard.add_pixel(pixel);
Material_hard_lin.add_pixel(pixel);
}
}
sys_lin.initialise();
sys_non.initialise();
Grad_t<dim> delF0;
delF0 << 1e-4, 5e-5, 5e-5, 0, 0, 0, 0, 0, 0;
constexpr Real cg_tol{1e-8}, newton_tol{1e-5};
constexpr Uint maxiter{CcoordOps::get_size(resolutions)*ipow(dim, secondOrder)*10};
constexpr bool verbose{false};
GradIncrements<dim> grads; grads.push_back(delF0);
DeprecatedSolverCG<dim> cg_lin{sys_lin, cg_tol, maxiter, bool(verbose)};
Eigen::ArrayXXd res2{deprecated_newton_cg(sys_lin, delF0, cg_lin, newton_tol, verbose).grad};
DeprecatedSolverCG<dim> cg_non{sys_non, cg_tol, maxiter, bool(verbose)};
Eigen::ArrayXXd res1{deprecated_newton_cg(sys_non, delF0, cg_non, newton_tol, verbose).grad};
BOOST_CHECK_LE(abs(res1-res2).mean(), cg_tol);
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
c4science · Help