Homec4science
Diffusion µSpectre (master) Browse

muSpectre/tests/split_test_intersection_sym_check.cc
f1f2d2107497master

Last Change
Raw File

split_test_intersection_sym_check.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 <fstream>
#include <boost/mpl/list.hpp>
#include <math.h>
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(symmetric_test);
BOOST_AUTO_TEST_CASE(checking_the_symmetry) {
constexpr Dim_t dim{twoD};
using Rcoord = Rcoord_t<dim>;
using Ccoord = Ccoord_t<dim>;
//typedef std::vector<Eigen::Array4d>::iterator Array_iter;
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 Real length {25};
constexpr Rcoord center{12.5, 12.5};
constexpr Rcoord halfside_neg{-3.1, -3.1};
constexpr Rcoord halfside_pos{+4.00, +4.00};
constexpr Real cg_tol{1e-8}, newton_tol{1e-5};
constexpr bool verbose{false};
std::vector<Rcoord> precipitate_vertices ;
precipitate_vertices.push_back ({center[0] + halfside_neg[0], center[1] + halfside_neg[1]});
precipitate_vertices.push_back ({center[0] + halfside_neg[0], center[1] + halfside_pos[1]});
precipitate_vertices.push_back ({center[0] + halfside_pos[0], center[1] + halfside_neg[1]});
precipitate_vertices.push_back ({center[0] + halfside_pos[0], center[1] + halfside_pos[1]});
Grad_t<dim> delF0_1;
delF0_1 << 0.002, 0.0,
0.001, 0.005;
////LOW RESOLUTION/////////////////////
constexpr int resolution_l{25};
constexpr Rcoord lengths_split_l{ length, length};
constexpr Ccoord resolutions_split_l{resolution_l, resolution_l};
auto fft_ptr_split_l{std::make_unique<FFTWEngine<dim>>(resolutions_split_l, ipow(dim, 2))};
auto proj_ptr_split_l{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr_split_l), lengths_split_l)};
CellSplit<dim, dim> sys_split_l(std::move(proj_ptr_split_l), SplittedCell::yes);
auto& Material_hard_split_l = Mat_t::make(sys_split_l, "hard", Young_hard, Poisson_hard);
auto& Material_soft_split_l = Mat_t::make(sys_split_l, "soft", Young_soft, Poisson_soft);
//assign precipitate materials:
RootNode<dim> precipitate_l (sys_split_l, precipitate_vertices);
//Extracting the intersected pixels and their correspondent intersection ratios:
auto precipitate_intersects_l = precipitate_l.get_intersected_pixels();
auto precipitate_intersection_ratios_l = precipitate_l.get_intersection_ratios();
Material_hard_split_l.add_split_pixels_precipitate(precipitate_intersects_l,
precipitate_intersection_ratios_l);
//std::vector<Real> assigned_ratio_l = sys_split_l.get_assigned_ratios();
auto incompleted_pixels_l = sys_split_l.make_incompleted_pixels();
for (auto && tup:incompleted_pixels_l){
auto && pix = std::get<0> (tup);
auto && ratio = std::get<1> (tup);
Material_soft_split_l.add_pixel_split(pix, ratio);
}
/* for (auto && tup: akantu::enumerate(sys_split_l)) {
auto && pixel = std::get<1>(tup);
auto iterator = std::get<0>(tup);
if (assigned_ratio_l[iterator] < 1.0){
Material_soft_split_l.add_pixel_split(pixel, 1.0-assigned_ratio_l[iterator]);
}
}*/
////HIGH RESOLUTION/////////////////////
constexpr int resolution_h{125};
constexpr Rcoord lengths_split_h{ length, length};
constexpr Ccoord resolutions_split_h{resolution_h, resolution_h};
auto fft_ptr_split_h{std::make_unique<FFTWEngine<dim>>(resolutions_split_h, ipow(dim, 2))};
auto proj_ptr_split_h{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr_split_h), lengths_split_h)};
CellSplit<dim, dim> sys_split_h(std::move(proj_ptr_split_h), SplittedCell::yes);
auto& Material_hard_split_h = Mat_t::make(sys_split_h, "hard", Young_hard, Poisson_hard);
auto& Material_soft_split_h = Mat_t::make(sys_split_h, "soft", Young_soft, Poisson_soft);
//assign precipitate materials:
RootNode<dim> precipitate_h (sys_split_h, precipitate_vertices);
//Extracting the intersected pixels and their correspondent intersection ratios:
auto precipitate_intersects_h= precipitate_h.get_intersected_pixels();
auto precipitate_intersection_ratios_h = precipitate_h.get_intersection_ratios();
Material_hard_split_h.add_split_pixels_precipitate(precipitate_intersects_h,
precipitate_intersection_ratios_h);
//std::vector<Real> assigned_ratio_h = sys_split_h.get_assigned_ratios();
/*
auto incompleted_pixels_h = sys_split_h.make_incompleted_pixels();
for (auto && tup:incompleted_pixels_h){
auto && pix = std::get<0> (tup);
auto && ratio = std::get<1> (tup);
Material_soft_split_h.add_pixel_split(pix, ratio);
}
*/
sys_split_h.complete_material_assignment(Material_soft_split_h);
/* for (auto && tup: akantu::enumerate(sys_split_h)) {
auto && pixel = std::get<1>(tup);
auto iterator = std::get<0>(tup);
if (assigned_ratio_h[iterator] < 1.0){
Material_soft_split_h.add_pixel_split(pixel, 1.0-assigned_ratio_h[iterator]);
}
}*/
//Finding equilbrium strain fileds with Sppectral method:
sys_split_l.initialise();
sys_split_h.initialise();
constexpr Uint maxiter_l{CcoordOps::get_size(resolutions_split_l)*ipow(dim, secondOrder)*10};
constexpr Uint maxiter_h{CcoordOps::get_size(resolutions_split_h)*ipow(dim, secondOrder)*10};
DeprecatedSolverCG<dim> cg_l{sys_split_l, cg_tol, maxiter_l, bool(verbose)};
auto output_l {deprecated_newton_cg(sys_split_l, delF0_1, cg_l, newton_tol, verbose)};
Eigen::ArrayXXd res_grad_l{output_l.grad};
Eigen::ArrayXXd res_P_l{output_l.stress};
DeprecatedSolverCG<dim> cg_h{sys_split_h, cg_tol, maxiter_h, bool(verbose)};
auto output_h {deprecated_newton_cg(sys_split_h, delF0_1, cg_h, newton_tol, verbose)};
Eigen::ArrayXXd res_grad_h{output_h.grad};
Eigen::ArrayXXd res_P_h{output_h.stress};
std::ofstream file_grad_l("grad_l.csv");
std::ofstream file_P_l("P_l.csv");
for (int j{0}; j < res_grad_l.cols(); j++){
file_grad_l << res_grad_l.col(j).transpose()<<"\n";
file_P_l << res_P_l.col(j).transpose()<<"\n";
}
std::ofstream file_grad_h("grad_h.csv");
std::ofstream file_P_h("P_h.csv");
for (int j{0}; j < res_grad_h.cols(); j++){
file_grad_h << res_grad_h.col(j).transpose()<<"\n";
file_P_h << res_P_h.col(j).transpose()<<"\n";
}
BOOST_CHECK_LE(abs(0),cg_tol);
}
BOOST_AUTO_TEST_SUITE_END();
}
c4science · Help