Homec4science
Diffusion µSpectre (intersection_octree) Browse

muSpectre/tests/test_patch_split_octree_cell_newton_cg_solver.cc
3b9c93d008c2intersection_octree

Last Change
Raw File

test_patch_split_octree_cell_newton_cg_solver.cc
View Options

/**
* @file test_patch_split_octree_cell_newton_cg_solver.cc
*
* @author Ali Faslfi <ali.faslfi@epfl.ch>
*
* @date 12 Jun 2018
*
* @brief Tests for split cells and octree material assignment
*
* Copyright © 2017 Till Ali Faslafi
*
* µ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 "solver/deprecated_solvers.hh"
#include "solver/deprecated_solver_cg.hh"
#include "solver/deprecated_solver_cg_eigen.hh"
#include "fft/fftw_engine.hh"
#include "fft/projection_finite_strain_fast.hh"
#include "materials/material_linear_elastic1.hh"
#include "common/iterators.hh"
#include "common/ccoord_operations.hh"
#include "common/common.hh"
#include "cell/cell_factory.hh"
#include "cell/cell_split.hh"
#include "common/intersection_octree.hh"
#include <boost/mpl/list.hpp>
#include <math.h>
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(split_cell_octree_tests);
BOOST_AUTO_TEST_CASE(octree_intersection) {
constexpr Dim_t dim{twoD};
using Rcoord = Rcoord_t<dim>;
using Ccoord = Ccoord_t<dim>;
//constexpr Dim_t dim{threeD};
using Mat_t = MaterialLinearElastic1<dim, dim>;
const Real contrast {10};
const Real Young_soft{1.0030648180242636}, Poisson_soft{0.29930675909878679};
const Real Young_hard{contrast * Young_soft}, Poisson_hard{0.29930675909878679};
constexpr Ccoord resolutions_split{11,11};
constexpr Rcoord lengths_split{1.1, 1.1};
auto fft_ptr_split_1{std::make_unique<FFTWEngine<dim>>(resolutions_split, ipow(dim, 2))};
auto proj_ptr_split_1{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr_split_1), lengths_split)};
auto fft_ptr_split_2{std::make_unique<FFTWEngine<dim>>(resolutions_split, ipow(dim, 2))};
auto proj_ptr_split_2{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr_split_2), lengths_split)};
CellSplit<dim, dim> sys_split_1(std::move(proj_ptr_split_1), SplittedCell::yes);
CellSplit<dim, dim> sys_split_2(std::move(proj_ptr_split_2), SplittedCell::yes);
auto& Material_hard_split_1 = Mat_t::make(sys_split_1, "hard", 10*Young_hard, Poisson_hard);
auto& Material_soft_split_1 = Mat_t::make(sys_split_1, "soft", Young_soft, Poisson_soft);
auto& Material_hard_split_2 = Mat_t::make(sys_split_2, "hard", 10*Young_hard, Poisson_hard);
auto& Material_soft_split_2 = Mat_t::make(sys_split_2, "soft", Young_soft, Poisson_soft);
//square precipitate vertices:
constexpr Rcoord center{0.55, 0.55};
constexpr Rcoord halfside{0.15, 0.15};
constexpr Real halfdiagon{sqrt( pow(halfside[0], 2) +
pow(halfside[1], 2) )};
//making the first precipitate
std::vector<Rcoord> precipitate_1_vertices ;
precipitate_1_vertices.push_back ({center[0] - halfside[0], center[1] - halfside[1]});
precipitate_1_vertices.push_back ({center[0] + halfside[0], center[1] - halfside[1]});
precipitate_1_vertices.push_back ({center[0] - halfside[0], center[1] + halfside[1]});
precipitate_1_vertices.push_back ({center[0] + halfside[0], center[1] + halfside[1]});
RootNode<dim> precipitate_1 (sys_split_1, precipitate_1_vertices);
auto precipitate_1_intersects = precipitate_1.get_intersected_pixels();
auto precipitate_1_intersection_ratios = precipitate_1.get_intersection_ratios();
//making the second precipitate
std::vector<Rcoord> precipitate_2_vertices ;
precipitate_2_vertices.push_back ({center[0] - halfdiagon, center[1] });
precipitate_2_vertices.push_back ({center[0] + halfdiagon, center[1] });
precipitate_2_vertices.push_back ({center[0] , center[1] - halfdiagon});
precipitate_2_vertices.push_back ({center[0] , center[1] + halfdiagon});
RootNode<dim> precipitate_2 (sys_split_2, precipitate_2_vertices);
auto precipitate_2_intersects = precipitate_2.get_intersected_pixels();
auto precipitate_2_intersection_ratios = precipitate_2.get_intersection_ratios();
//assign precipitate materials:
typedef std::vector<Ccoord>::iterator Ccoord_iter_1, Ccoord_iter_2;
typedef std::vector<Real>::iterator Real_iter_1, Real_iter_2;
for (std::pair<Ccoord_iter_1, Real_iter_1> it(precipitate_1_intersects.begin(), precipitate_1_intersection_ratios.begin()) ;
it.first != precipitate_1_intersects.end(); ++it.first , ++it.second ){
Material_hard_split_1.add_pixel_split(*it.first, *it.second);
Material_soft_split_1.add_pixel_split(*it.first, 1.0-*it.second);
}
for (std::pair<Ccoord_iter_2, Real_iter_2> it(precipitate_2_intersects.begin(), precipitate_2_intersection_ratios.begin()) ;
it.first != precipitate_2_intersects.end(); ++it.first , ++it.second ){
Material_hard_split_2.add_pixel_split(*it.first, *it.second);
Material_soft_split_2.add_pixel_split(*it.first, 1.0-*it.second);
}
//complete material assignment
//sys_split_1.complete_material_assignment<Mat_t>(Material_soft_split);
std::vector<Real> assigned_ratio_1 = sys_split_1.get_intersected_ratios();
int iterator = 0 ;
for (auto && tup: akantu::enumerate(sys_split_1)) {
auto && pixel = std::get<1>(tup);
if (assigned_ratio_1[iterator] != 1.0){
Material_soft_split_1.add_pixel_split(pixel, 1.0-assigned_ratio_1[iterator]);
iterator++;
}
}
std::vector<Real> assigned_ratio_2 = sys_split_2.get_intersected_ratios();
iterator = 0 ;
for (auto && tup: akantu::enumerate(sys_split_2)) {
auto && pixel = std::get<1>(tup);
if (assigned_ratio_2[iterator] != 1.0){
Material_soft_split_2.add_pixel_split(pixel, 1.0-assigned_ratio_2[iterator]);
iterator++;
}
}
sys_split_1.initialise();
sys_split_2.initialise();
Eigen::Rotation2D<Real> rotation_1_2( pi/4);
Eigen::Rotation2D<Real> rotation_2_1(-pi/4);
Grad_t<dim> delF0_1, delF0_2;
delF0_1 << 0, 1.0, 0, 0;
delF0_2 = rotation_1_2 * delF0_1;
constexpr Real cg_tol_1{1e-8}, newton_tol_1{1e-5};
constexpr Uint maxiter_1{CcoordOps::get_size(resolutions_split)*ipow(dim, secondOrder)*10};
constexpr bool verbose_1{false};
GradIncrements<dim> grads_1; grads_1.push_back(delF0_1);
DeprecatedSolverCG<dim> cg_1{sys_split_1, cg_tol_1, maxiter_1, bool(verbose_1)};
Eigen::ArrayXXd res1{deprecated_de_geus(sys_split_1, grads_1, cg_1, newton_tol_1, verbose_1)[0].grad};
DeprecatedSolverCGEigen<dim> cg_2{sys_split_2, cg_tol_1, maxiter_1, bool(verbose_1)};
Eigen::ArrayXXd res2{deprecated_de_geus(sys_split_2, grads_1, cg_2, newton_tol_1, verbose_1)[0].grad};
BOOST_CHECK_LE(abs(0), cg_tol_1);
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
c4science · Help