Homec4science
Diffusion µSpectre (feat/laminate_pixels) Browse

muSpectre/tests/split_test_intersection_error_induced.cc
5e584d008633feat/laminate_pixels

Last Change
Raw File

split_test_intersection_error_induced.cc
View Options

/**
* @file test_intersection_error_induced.cc
*
* @author Ali Faslfi <ali.faslfi@epfl.ch>
*
* @date 21 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/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_error_test);
BOOST_AUTO_TEST_CASE(area_calculation_and_eq_solution_of_55vs11_pixels) {
constexpr Dim_t dim{twoD};
using Rcoord = Rcoord_t<dim>;
using Ccoord = Ccoord_t<dim>;
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};
const Real machine_precision{1e-13};
constexpr Real length {11};
constexpr int high_res{55},low_res{11};
constexpr Rcoord lengths_split{length, length};
constexpr Ccoord resolutions_split_high_res{high_res, high_res};
constexpr Ccoord resolutions_split_low_res{low_res, low_res};
auto fft_ptr_split_high_res{std::make_unique<FFTWEngine<dim>>(resolutions_split_high_res, ipow(dim, 2))};
auto fft_ptr_split_low_res {std::make_unique<FFTWEngine<dim>>(resolutions_split_low_res, ipow(dim, 2))};
auto proj_ptr_split_high_res{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>
(std::move(fft_ptr_split_high_res), lengths_split)};
auto proj_ptr_split_low_res{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>
(std::move(fft_ptr_split_low_res), lengths_split)};
CellSplit<dim, dim> sys_split_high_res(std::move(proj_ptr_split_high_res), SplittedCell::yes);
CellSplit<dim, dim> sys_split_low_res(std::move(proj_ptr_split_low_res), SplittedCell::yes);
auto& Material_hard_split_high_res = Mat_t::make(sys_split_high_res, "hard", Young_hard, Poisson_hard);
auto& Material_soft_split_high_res = Mat_t::make(sys_split_high_res, "soft", Young_soft, Poisson_soft);
auto& Material_hard_split_low_res = Mat_t::make(sys_split_low_res, "hard", Young_hard, Poisson_hard);
auto& Material_soft_split_low_res = Mat_t::make(sys_split_low_res, "soft", Young_soft, Poisson_soft);
constexpr Rcoord center{5.5, 5.5};
constexpr Rcoord halfside1{2.3, 2.3};
constexpr Rcoord halfside2{2.4, 2.4};
std::vector<Rcoord> precipitate_vertices ;
precipitate_vertices.push_back ({center[0] - halfside1[0], center[1] - halfside1[1]});
precipitate_vertices.push_back ({center[0] + halfside2[0], center[1] - halfside1[1]});
precipitate_vertices.push_back ({center[0] - halfside1[0], center[1] + halfside2[1]});
precipitate_vertices.push_back ({center[0] + halfside2[0], center[1] + halfside2[1]});
//analyzing the intersection of the preicipitate with the pixels
RootNode<dim> precipitate_high_res (sys_split_high_res, precipitate_vertices);
//Extracting the intersected pixels and their correspondent intersection ratios:
auto precipitate_high_res_intersects = precipitate_high_res.get_intersected_pixels();
auto precipitate_high_res_intersection_ratios = precipitate_high_res.get_intersection_ratios();
//assign material to the pixels which have intersection with the precipitate
for (auto tup:akantu::zip(precipitate_high_res_intersects, precipitate_high_res_intersection_ratios)){
auto pix = std::get<0> (tup);
auto ratio = std::get<1>(tup);
Material_hard_split_high_res.add_pixel_split(pix,ratio);
Material_soft_split_high_res.add_pixel_split(pix, 1.0-ratio);
}
//assign material to the rest of the pixels:
std::vector<Real> assigned_ratio_high_res = sys_split_high_res.get_assigned_ratios();
for (auto && tup: akantu::enumerate(sys_split_high_res)) {
auto && pixel = std::get<1>(tup);
auto iterator = std::get<0>(tup);
if (assigned_ratio_high_res[iterator] < 1.0){
Material_soft_split_high_res.add_pixel_split(pixel, 1.0-assigned_ratio_high_res[iterator]);
}
}
sys_split_high_res.initialise();
//analyzing the intersection of the preicipitate with the pixels
RootNode<dim> precipitate_low_res (sys_split_low_res, precipitate_vertices);
//Extracting the intersected pixels and their correspondent intersection ratios:
auto precipitate_low_res_intersects = precipitate_low_res.get_intersected_pixels();
auto precipitate_low_res_intersection_ratios = precipitate_low_res.get_intersection_ratios();
//assign material to the pixels which have intersection with the precipitate
for (auto tup:akantu::zip(precipitate_low_res_intersects, precipitate_low_res_intersection_ratios)){
auto pix = std::get<0> (tup);
auto ratio = std::get<1>(tup);
Material_hard_split_low_res.add_pixel_split(pix, ratio);
Material_soft_split_low_res.add_pixel_split(pix, 1.0-ratio);
}
//assign material to the rest of the pixels:
std::vector<Real> assigned_ratio_low_res = sys_split_low_res.get_assigned_ratios();
for (auto && tup: akantu::enumerate(sys_split_low_res)) {
auto && pixel = std::get<1>(tup);
auto iterator = std::get<0>(tup);
if (assigned_ratio_low_res[iterator] < 1.0){
Material_soft_split_low_res.add_pixel_split(pixel, 1.0-assigned_ratio_low_res[iterator]);
}
}
sys_split_low_res.initialise();
//Calculating the area of the precepetiate from the intersected pixels:
Real area_preticipitate_low_res = 0.0 ,area_preticipitate_high_res = 0.0 ;
for (auto&& precepetiate_area_low_res: precipitate_low_res_intersection_ratios){
area_preticipitate_low_res += precepetiate_area_low_res * sys_split_low_res.get_pixel_volume();
}
for (auto&& precepetiate_area_high_res: precipitate_high_res_intersection_ratios){
area_preticipitate_high_res += precepetiate_area_high_res * sys_split_high_res.get_pixel_volume();
}
BOOST_CHECK_LE((abs(area_preticipitate_high_res -
area_preticipitate_low_res)
/area_preticipitate_low_res),
machine_precision);
/////////////////////////////////////////////////////////////////////////////////////////////////
/*
//Finding equilbrium strain fileds with Sppectral method:
sys_split_low_res.initialise();
sys_split_high_res.initialise();
Grad_t<dim> delF0_1;
delF0_1 << 1.0, 0.0,
0.0, 1.0;
constexpr Real cg_tol_low_res{1e-8}, newton_tol_low_res{1e-5};
constexpr Uint maxiter_low_res{CcoordOps::get_size(resolutions_split_low_res)*ipow(dim, secondOrder)*10};
constexpr bool verbose_low_res{false};
constexpr Real cg_tol_high_res{1e-8}, newton_tol_high_res{1e-5};
constexpr Uint maxiter_high_res{CcoordOps::get_size(resolutions_split_high_res)*ipow(dim, secondOrder)*10};
constexpr bool verbose_high_res{false};
DeprecatedSolverCG<dim> cg_low_res{sys_split_low_res, cg_tol_low_res, maxiter_low_res, bool(verbose_low_res)};
Eigen::ArrayXXd res_low_res{deprecated_newton_cg(sys_split_low_res, delF0_1, cg_low_res, newton_tol_low_res, verbose_low_res).grad};
DeprecatedSolverCGEigen<dim> cg_high_res{sys_split_high_res, cg_tol_high_res, maxiter_high_res, bool(verbose_high_res)};
Eigen::ArrayXXd res_high_res{deprecated_newton_cg(sys_split_high_res, delF0_1, cg_high_res, newton_tol_high_res, verbose_high_res).grad};
/std::vector<std::vector<Eigen::Array4d>> result_low_res{}, result_high_res_filtered{}, result_comparison{} ;
std::vector<Eigen::Array4d> temp_vec_low_res{}, temp_vec_high_res{}, temp_vec_result_comparison{};
for (int i{0} ; i < high_res ; i++){
if (i%5 == 0){
for(int j{0} ; j < high_res ; j++){
if(j%5 == 0){
temp_vec_high_res.push_back(res_high_res.col(i* high_res +j));
}
}
result_high_res_filtered.push_back(temp_vec_high_res);
temp_vec_high_res.clear();
}
}
for (int i{0} ; i < low_res ; i++){
for(int j{0} ; j < low_res ; j++){
temp_vec_low_res.push_back(res_low_res.col(i* low_res +j));
}
result_low_res.push_back(temp_vec_low_res);
temp_vec_low_res.clear();
}
for (unsigned int i{0} ; i < result_low_res.size(); i++){
for (unsigned int j{0} ; j < result_low_res.size(); j++){
temp_vec_result_comparison.push_back(result_high_res_filtered[i][j]-result_low_res[i][j]);
}
result_comparison.push_back(temp_vec_result_comparison);
temp_vec_result_comparison.clear();
}*
BOOST_CHECK_LE(abs(0),cg_tol_low_res);*/
}
BOOST_AUTO_TEST_SUITE_END();
}
c4science · Help