diff --git a/src/common/geometry.hh b/src/common/geometry.hh index cfb81d8..72c58ff 100644 --- a/src/common/geometry.hh +++ b/src/common/geometry.hh @@ -1,507 +1,507 @@ /** * @file geometry.hh * * @author Till Junge * * @date 18 Apr 2018 * * @brief Geometric calculation helpers * * Copyright © 2018 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 "common/common.hh" #include "common/tensor_algebra.hh" #include "common/eigen_tools.hh" #include #include #include namespace muSpectre { /** * The rotation matrices depend on the order in which we rotate * around different axes. See [[ * https://en.wikipedia.org/wiki/Euler_angles#Rotation_matrix ]] to * find the matrices */ enum class RotationOrder {Z, XZXEuler, XYXEuler, YXYEuler, YZYEuler, ZYZEuler, ZXZEuler, XZYTaitBryan, XYZTaitBryan, YXZTaitBryan, YZXTaitBryan, ZYXTaitBryan, ZXYTaitBryan}; namespace internal { template struct DefaultOrder { constexpr static RotationOrder value{RotationOrder::ZXYTaitBryan}; }; template <> struct DefaultOrder { constexpr static RotationOrder value{RotationOrder::Z}; }; } // internal template class RotatorBase{ public: using RotMat_t = Eigen::Matrix; //! Default constructor RotatorBase() = delete; RotatorBase(const RotMat_t rotation_matrix_input): rot_mat{rotation_matrix_input} {} //! Copy constructor RotatorBase(const RotatorBase &other) = default; //! Move constructor RotatorBase(RotatorBase &&other) = default; //! Destructor virtual ~RotatorBase() = default; //! Copy assignment operator RotatorBase& operator=(const RotatorBase &other) = default; //! Move assignment operator RotatorBase& operator=(RotatorBase &&other) = default; /** * Applies the rotation into the frame define my the rotation * matrix * * @param input is a first-, second-, or fourth-rank tensor * (column vector, square matrix, or T4Matrix, or a Eigen::Map of * either of these, or an expression that evaluates into any of * these) */ template inline decltype(auto) rotate(In_t && input); /** * Applies the rotation back out from the frame define my the * rotation matrix * * @param input is a first-, second-, or fourth-rank tensor * (column vector, square matrix, or T4Matrix, or a Eigen::Map of * either of these, or an expression that evaluates into any of * these) */ template inline decltype(auto) rotate_back(In_t && input); const RotMat_t & get_rot_mat() const {return this->rot_mat;} protected: EIGEN_MAKE_ALIGNED_OPERATOR_NEW; RotMat_t rot_mat; }; template ::value> class RotatorAngle: public RotatorBase { static_assert(((Dim == twoD) and (Order == RotationOrder::Z)) or ((Dim == threeD) and (Order != RotationOrder::Z)), "In 2d, only order 'Z' makes sense. In 3d, it doesn't"); public: using Parent = RotatorBase; using Angles_t = Eigen::Matrix; using RotMat_t = Eigen::Matrix; //! Default constructor RotatorAngle() = delete; RotatorAngle(const Eigen::Ref & angles_inp) :Parent(this->compute_rotation_matrix_angle(angles_inp)), angles{angles_inp} {} //! Copy constructor RotatorAngle(const RotatorAngle &other) = default; //! Move constructor RotatorAngle(RotatorAngle &&other) = default; //! Destructor virtual ~RotatorAngle() = default; //! Copy assignment operator RotatorAngle& operator=(const RotatorAngle &other) = default; //! Move assignment operator RotatorAngle& operator=(RotatorAngle &&other) = default; protected: inline RotMat_t compute_rotation_matrix_angle(Angles_t angles); inline RotMat_t compute_this_rotation_matrix_angle(){ return compute_rotation_matrix(this->angles); }; EIGEN_MAKE_ALIGNED_OPERATOR_NEW; Angles_t angles; private: }; /* ---------------------------------------------------------------------- */ namespace internal { template struct RotationMatrixComputerAngle {}; template struct RotationMatrixComputerAngle { constexpr static Dim_t Dim{twoD}; using RotMat_t = typename RotatorAngle::RotMat_t; using Angles_t = typename RotatorAngle::Angles_t; inline static decltype(auto) compute(const Eigen::Ref & angles) { static_assert(Order == RotationOrder::Z, "Two-d rotations can only be around the z axis"); return RotMat_t(Eigen::Rotation2Dd(angles(0))); } }; template struct RotationMatrixComputerAngle { constexpr static Dim_t Dim{threeD}; using RotMat_t = typename RotatorAngle::RotMat_t; using Angles_t = typename RotatorAngle::Angles_t; inline static decltype(auto) compute(const Eigen::Ref & angles) { static_assert(Order != RotationOrder::Z, "three-d rotations cannot only be around the z axis"); switch (Order) { case RotationOrder::ZXZEuler: { return RotMat_t((Eigen::AngleAxisd(angles(0), Eigen::Vector3d::UnitZ()) * Eigen::AngleAxisd(angles(1), Eigen::Vector3d::UnitX()) * Eigen::AngleAxisd(angles(2), Eigen::Vector3d::UnitZ()))); break; } case RotationOrder::ZXYTaitBryan: { return RotMat_t((Eigen::AngleAxisd(angles(0), Eigen::Vector3d::UnitZ()) * Eigen::AngleAxisd(angles(1), Eigen::Vector3d::UnitX()) * Eigen::AngleAxisd(angles(2), Eigen::Vector3d::UnitY()))); } default: { throw std::runtime_error("not yet implemented."); } } } }; } // internal /* ---------------------------------------------------------------------- */ template auto RotatorAngle:: compute_rotation_matrix_angle(Angles_t angles) -> RotMat_t { return internal::RotationMatrixComputerAngle::compute(angles); } /* ---------------------------------------------------------------------- */ /** * this class is used to make the vector a aligned to the vec b by means of a rotation system, the input for the constructor is the vector itself and the functions rotate and rotate back would be available as they exist in the parent class (RotatorBase) nad can be used in order to do the functionality of the class */ template class RotatorTwoVec: public RotatorBase { public: using Parent = RotatorBase; using Vec_t = Eigen::Matrix; using RotMat_t = Eigen::Matrix; //! Default constructor RotatorTwoVec() = delete; RotatorTwoVec(Vec_t vec_a_inp, Vec_t vec_b_inp) :Parent(this->compute_rotation_matrix_TwoVec(vec_a_inp, vec_b_inp)), vec_ref{vec_a_inp}, vec_des{vec_b_inp} {} //! Copy constructor RotatorTwoVec(const RotatorTwoVec &other) = default; //! Move constructor RotatorTwoVec(RotatorTwoVec &&other) = default; //! Destructor virtual ~RotatorTwoVec() = default; //! Copy assignment operator RotatorTwoVec& operator=(const RotatorTwoVec &other) = default; //! Move assignment operator RotatorTwoVec& operator=(RotatorTwoVec &&other) = default; protected: inline RotMat_t compute_rotation_matrix_TwoVec(Vec_t vec_ref, Vec_t vec_des); inline RotMat_t compute_this_rotation_matrix_TwoVec(){ return compute_rotation_matrix(this->vec_ref, this->vec_des); } EIGEN_MAKE_ALIGNED_OPERATOR_NEW; Vec_t vec_ref, vec_des; private: }; /* ---------------------------------------------------------------------- */ namespace internal { template struct RotationMatrixComputerTwoVec {}; template<> struct RotationMatrixComputerTwoVec { constexpr static Dim_t Dim{twoD}; using RotMat_t = typename RotatorTwoVec::RotMat_t; using Vec_t = typename RotatorTwoVec::Vec_t; inline static RotMat_t compute(Vec_t vec_ref, Vec_t vec_des) { Real v_ref_norm = sqrt(vec_ref(0)*vec_ref(0) + vec_ref(1)*vec_ref(1)); Real v_des_norm = sqrt(vec_des(0)*vec_des(0) + vec_des(1)*vec_des(1)); RotMat_t ret_mat ; ret_mat(0,0) = ret_mat(1,1) = (((vec_ref(0)/v_ref_norm)* (vec_des(0)/v_des_norm))+ ((vec_des(1)/v_des_norm)* (vec_ref(1)/v_ref_norm))); ret_mat(1,0) = (((vec_ref(0)/v_ref_norm)* (vec_des(1)/v_des_norm))- ((vec_des(0)/v_des_norm)* (vec_ref(1)/v_ref_norm))); ret_mat(0,1 ) = -ret_mat(1,0); return ret_mat; } }; template<> struct RotationMatrixComputerTwoVec { constexpr static Dim_t Dim{threeD}; using RotMat_t = typename RotatorTwoVec::RotMat_t; using Vec_t = typename RotatorTwoVec::Vec_t; inline static RotMat_t compute(Vec_t vec_ref, Vec_t vec_des) { return Eigen::Quaternion::FromTwoVectors(vec_ref, vec_des). normalized().toRotationMatrix(); } }; } //internal /* ---------------------------------------------------------------------- */ template auto RotatorTwoVec:: compute_rotation_matrix_TwoVec(Vec_t vec_ref, Vec_t vec_des) -> RotMat_t { return internal::RotationMatrixComputerTwoVec::compute(vec_ref, vec_des); } /* ---------------------------------------------------------------------- */ /** * this class is used to make a vector aligned to x-axis of the coordinate system, the input for the constructor is the vector itself and the functions rotate and rotate back would be available as they exist in the parent class (RotatorBase) nad can be used in order to do the functionality of the class */ template class RotatorNormal: public RotatorBase { public: using Parent = RotatorBase; - using Vec_t = Eigen::Matrix; + using Vec_t = Eigen::Matrix; using RotMat_t = Eigen::Matrix; //! Default constructor RotatorNormal() = delete; RotatorNormal(Vec_t vec) :Parent(this->compute_rotation_matrix_normal(vec)), vec{vec} {} //! Copy constructor RotatorNormal(const RotatorNormal &other) = default; //! Move constructor RotatorNormal(RotatorNormal &&other) = default; //! Destructor virtual ~RotatorNormal() = default; //! Copy assignment operator RotatorNormal& operator=(const RotatorNormal &other) = default; //! Move assignment operator RotatorNormal& operator=(RotatorNormal &&other) = default; protected: inline RotMat_t compute_rotation_matrix_normal(Vec_t vec); inline RotMat_t compute_this_rotation_matrix_normal(){ return compute_rotation_matrix_normal(this->vec); } EIGEN_MAKE_ALIGNED_OPERATOR_NEW; Vec_t vec; private: }; /* ---------------------------------------------------------------------- */ namespace internal { template struct RotationMatrixComputerNormal {}; template<> struct RotationMatrixComputerNormal { constexpr static Dim_t Dim{twoD}; using RotMat_t = typename RotatorTwoVec::RotMat_t; using Vec_t = typename RotatorTwoVec::Vec_t; inline static RotMat_t compute(Vec_t vec) { Real vec_norm = sqrt(vec(0)*vec(0) + vec(1)*vec(1)); Vec_t x; x<<1.0, 0.0; RotMat_t ret_mat; ret_mat(0,0) = ret_mat(1,1) = ((vec(0)/vec_norm)* x(0)); ret_mat(1,0) = -(-(vec(1)/vec_norm)* x(0)); ret_mat(0,1) = - ret_mat(1,0); return ret_mat; } }; template<> struct RotationMatrixComputerNormal { constexpr static Dim_t Dim{threeD}; using RotMat_t = typename RotatorTwoVec::RotMat_t; using Vec_t = typename RotatorTwoVec::Vec_t; inline static RotMat_t compute(Vec_t vec) { Real eps = 0.1 ; Vec_t vec1 = vec/vec.norm(); Vec_t x(Vec_t::UnitX()); Vec_t y(Vec_t::UnitY()); Vec_t n_x = vec1.cross(x); Vec_t vec2 = ((n_x.norm()>eps) * n_x + (1-(n_x.norm()>eps))*(vec1.cross(y))); Vec_t vec3 = vec1.cross(vec2); RotMat_t ret_mat; ret_mat < auto RotatorNormal:: compute_rotation_matrix_normal(Vec_t vec) -> RotMat_t{ return internal::RotationMatrixComputerNormal::compute(vec); } /* ---------------------------------------------------------------------- */ namespace internal { template struct RotationHelper { }; /* ---------------------------------------------------------------------- */ template <> struct RotationHelper { template inline static decltype(auto) rotate(In_t && input, Rot_t && R) { return R * input; } }; /* ---------------------------------------------------------------------- */ template <> struct RotationHelper { template inline static decltype(auto) rotate(In_t && input, Rot_t && R) { return R * input * R.transpose(); } }; /* ---------------------------------------------------------------------- */ template <> struct RotationHelper { template inline static decltype(auto) rotate(In_t && input, Rot_t && R) { constexpr Dim_t Dim{EigenCheck::tensor_dim::value}; auto && rotator_forward{Matrices::outer_under(R.transpose(), R.transpose())}; auto && rotator_back = Matrices::outer_under(R, R); // unclear behaviour. When I return this value as an // expression, clange segfaults or returns an uninitialised // tensor return T4Mat(rotator_back * input * rotator_forward); } }; } /* ---------------------------------------------------------------------- */ template template auto RotatorBase::rotate(In_t && input) -> decltype(auto) { constexpr Dim_t tensor_rank{EigenCheck::tensor_rank::value}; return internal::RotationHelper:: rotate(std::forward(input), this->rot_mat); } /* ---------------------------------------------------------------------- */ template template auto RotatorBase::rotate_back(In_t && input) -> decltype(auto) { constexpr Dim_t tensor_rank{EigenCheck::tensor_rank::value}; return internal::RotationHelper:: rotate(std::forward(input), this->rot_mat.transpose()); } } // muSpectre diff --git a/tests/asplit_test_intersection_error_induced.cc b/tests/asplit_test_intersection_error_induced.cc new file mode 100644 index 0000000..1b907c6 --- /dev/null +++ b/tests/asplit_test_intersection_error_induced.cc @@ -0,0 +1,215 @@ +/** + * @file test_intersection_error_induced.cc + * + * @author Ali Faslfi + * + * @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 +#include + +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; + using Ccoord = Ccoord_t; + + using Mat_t = MaterialLinearElastic1; + 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>(resolutions_split_high_res, ipow(dim, 2))}; + auto fft_ptr_split_low_res {std::make_unique>(resolutions_split_low_res, ipow(dim, 2))}; + + auto proj_ptr_split_high_res{std::make_unique> + (std::move(fft_ptr_split_high_res), lengths_split)}; + auto proj_ptr_split_low_res{std::make_unique> + (std::move(fft_ptr_split_low_res), lengths_split)}; + + CellSplit sys_split_high_res(std::move(proj_ptr_split_high_res), SplitCell::simple); + CellSplit sys_split_low_res(std::move(proj_ptr_split_low_res), SplitCell::simple); + + + 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 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 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 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 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 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 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 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 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> result_low_res{}, result_high_res_filtered{}, result_comparison{} ; + std::vector 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(); +} diff --git a/tests/asplit_test_intersection_sym_check.cc b/tests/asplit_test_intersection_sym_check.cc new file mode 100644 index 0000000..b168974 --- /dev/null +++ b/tests/asplit_test_intersection_sym_check.cc @@ -0,0 +1,191 @@ +/** + * @file test_intersection_error_induced.cc + * + * @author Ali Faslfi + * + * @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 +#include +#include + +namespace muSpectre { + + BOOST_AUTO_TEST_SUITE(symmetric_test); + + BOOST_AUTO_TEST_CASE(checking_the_symmetry) { + constexpr Dim_t dim{twoD}; + using Rcoord = Rcoord_t; + using Ccoord = Ccoord_t; + //typedef std::vector::iterator Array_iter; + using Mat_t = MaterialLinearElastic1; + 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 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 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>(resolutions_split_l, ipow(dim, 2))}; + auto proj_ptr_split_l{std::make_unique>(std::move(fft_ptr_split_l), lengths_split_l)}; + CellSplit sys_split_l(std::move(proj_ptr_split_l), SplitCell::simple); + + 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); + sys_split_l.make_automatic_precipitate(precipitate_vertices, Material_hard_split_l); + sys_split_l.complete_material_assignment(Material_soft_split_l); + //assign precipitate materials: + /* RootNode 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 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>(resolutions_split_h, ipow(dim, 2))}; + auto proj_ptr_split_h{std::make_unique>(std::move(fft_ptr_split_h), lengths_split_h)}; + CellSplit sys_split_h(std::move(proj_ptr_split_h), SplitCell::simple); + + + 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); + + sys_split_h.make_automatic_precipitate(precipitate_vertices, Material_hard_split_h); + sys_split_h.complete_material_assignment(Material_soft_split_h); + //assign precipitate materials: + /*RootNode 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 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); + }*/ + + /* 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 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 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(); +}