diff --git a/bin/demonstrator1.cc b/bin/demonstrator1.cc index 504d4c3..4810549 100644 --- a/bin/demonstrator1.cc +++ b/bin/demonstrator1.cc @@ -1,142 +1,144 @@ /** * @file demonstrator1.cc * * @author Till Junge * * @date 03 Jan 2018 * * @brief larger problem to show off * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include #include #include "external/cxxopts.hpp" #include "common/muSpectre_common.hh" #include #include "cell/cell_factory.hh" #include "materials/material_linear_elastic1.hh" #include "solver/deprecated_solvers.hh" #include "solver/deprecated_solver_cg.hh" using opt_ptr = std::unique_ptr; opt_ptr parse_args(int argc, char ** argv) { opt_ptr options = std::make_unique(argv[0], "Tests MPI fft scalability"); try { options->add_options()("0,N0", "number of rows", cxxopts::value(), "N0")("h,help", "print help")( "positional", "Positional arguments: these are the arguments that are entered " "without an option", cxxopts::value>()); options->parse_positional(std::vector{"N0", "positional"}); options->parse(argc, argv); if (options->count("help")) { std::cout << options->help({"", "Group"}) << std::endl; exit(0); } if (options->count("N0") != 1) { throw cxxopts::OptionException("Parameter N0 missing"); } else if ((*options)["N0"].as() % 2 != 1) { throw cxxopts::OptionException("N0 must be odd"); } else if (options->count("positional") > 0) { throw cxxopts::OptionException("There are too many positional arguments"); } } catch (const cxxopts::OptionException & e) { std::cout << "Error parsing options: " << e.what() << std::endl; exit(1); } return options; } +using namespace muFFT; +using namespace muGrid; using namespace muSpectre; int main(int argc, char * argv[]) { banner("demonstrator1", 2018, "Till Junge "); auto options{parse_args(argc, argv)}; auto & opt{*options}; const Dim_t size{opt["N0"].as()}; constexpr Real fsize{1.}; constexpr Dim_t dim{3}; const Dim_t nb_dofs{ipow(size, dim) * ipow(dim, 2)}; std::cout << "Number of dofs: " << nb_dofs << std::endl; constexpr Formulation form{Formulation::finite_strain}; const Rcoord_t lengths{CcoordOps::get_cube(fsize)}; const Ccoord_t resolutions{CcoordOps::get_cube(size)}; auto cell{make_cell(resolutions, lengths, form)}; constexpr Real E{1.0030648180242636}; constexpr Real nu{0.29930675909878679}; using Material_t = MaterialLinearElastic1; auto & Material_soft{Material_t::make(cell, "soft", E, nu)}; auto & Material_hard{Material_t::make(cell, "hard", 10 * E, nu)}; int counter{0}; for (const auto && pixel : cell) { int sum = 0; for (Dim_t i = 0; i < dim; ++i) { sum += pixel[i] * 2 / resolutions[i]; } if (sum == 0) { Material_hard.add_pixel(pixel); counter++; } else { Material_soft.add_pixel(pixel); } } std::cout << counter << " Pixel out of " << cell.size() << " are in the hard material" << std::endl; cell.initialise(FFT_PlanFlags::measure); constexpr Real newton_tol{1e-4}; constexpr Real cg_tol{1e-7}; const Uint maxiter = nb_dofs; Grad_t DeltaF{Grad_t::Zero()}; DeltaF(0, 1) = .1; Dim_t verbose{1}; auto start = std::chrono::high_resolution_clock::now(); GradIncrements grads{DeltaF}; DeprecatedSolverCG cg{cell, cg_tol, maxiter, static_cast(verbose)}; deprecated_newton_cg(cell, grads, cg, newton_tol, verbose); std::chrono::duration dur = std::chrono::high_resolution_clock::now() - start; std::cout << "Resolution time = " << dur.count() << "s" << std::endl; return 0; } diff --git a/bin/demonstrator_dynamic_solve.cc b/bin/demonstrator_dynamic_solve.cc index feb394d..d19defd 100644 --- a/bin/demonstrator_dynamic_solve.cc +++ b/bin/demonstrator_dynamic_solve.cc @@ -1,142 +1,145 @@ /** * @file demonstrator1.cc * * @author Till Junge * * @date 03 Jan 2018 * * @brief larger problem to show off * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include #include #include "external/cxxopts.hpp" #include "common/muSpectre_common.hh" -#include #include "cell/cell_factory.hh" #include "materials/material_linear_elastic1.hh" #include "solver/solvers.hh" #include "solver/solver_cg.hh" +#include + using opt_ptr = std::unique_ptr; opt_ptr parse_args(int argc, char ** argv) { opt_ptr options = std::make_unique(argv[0], "Tests MPI fft scalability"); try { options->add_options()("0,N0", "number of rows", cxxopts::value(), "N0")("h,help", "print help")( "positional", "Positional arguments: these are the arguments that are entered " "without an option", cxxopts::value>()); options->parse_positional(std::vector{"N0", "positional"}); options->parse(argc, argv); if (options->count("help")) { std::cout << options->help({"", "Group"}) << std::endl; exit(0); } if (options->count("N0") != 1) { throw cxxopts::OptionException("Parameter N0 missing"); } else if ((*options)["N0"].as() % 2 != 1) { throw cxxopts::OptionException("N0 must be odd"); } else if (options->count("positional") > 0) { throw cxxopts::OptionException("There are too many positional arguments"); } } catch (const cxxopts::OptionException & e) { std::cout << "Error parsing options: " << e.what() << std::endl; exit(1); } return options; } using namespace muSpectre; +using namespace muGrid; +using namespace muFFT; int main(int argc, char * argv[]) { banner("demonstrator1", 2018, "Till Junge "); auto options{parse_args(argc, argv)}; auto & opt{*options}; const Dim_t size{opt["N0"].as()}; constexpr Real fsize{1.}; constexpr Dim_t dim{3}; const Dim_t nb_dofs{ipow(size, dim) * ipow(dim, 2)}; std::cout << "Number of dofs: " << nb_dofs << std::endl; constexpr Formulation form{Formulation::finite_strain}; const Rcoord_t lengths{CcoordOps::get_cube(fsize)}; const Ccoord_t resolutions{CcoordOps::get_cube(size)}; auto cell{make_cell(resolutions, lengths, form)}; constexpr Real E{1.0030648180242636}; constexpr Real nu{0.29930675909878679}; using Material_t = MaterialLinearElastic1; auto & Material_soft{Material_t::make(cell, "soft", E, nu)}; auto & Material_hard{Material_t::make(cell, "hard", 10 * E, nu)}; int counter{0}; for (const auto && pixel : cell) { int sum = 0; for (Dim_t i = 0; i < dim; ++i) { sum += pixel[i] * 2 / resolutions[i]; } if (sum == 0) { Material_hard.add_pixel(pixel); counter++; } else { Material_soft.add_pixel(pixel); } } std::cout << counter << " Pixel out of " << cell.size() << " are in the hard material" << std::endl; cell.initialise(FFT_PlanFlags::measure); constexpr Real newton_tol{1e-4}; constexpr Real cg_tol{1e-7}; const Uint maxiter = nb_dofs; Eigen::MatrixXd DeltaF{Eigen::MatrixXd::Zero(dim, dim)}; DeltaF(0, 1) = .1; Dim_t verbose{1}; auto start = std::chrono::high_resolution_clock::now(); LoadSteps_t loads{DeltaF}; SolverCG cg{cell, cg_tol, maxiter, static_cast(verbose)}; newton_cg(cell, loads, cg, newton_tol, verbose); std::chrono::duration dur = std::chrono::high_resolution_clock::now() - start; std::cout << "Resolution time = " << dur.count() << "s" << std::endl; return 0; } diff --git a/bin/small_case.cc b/bin/small_case.cc index 1a6e8c3..485ebfd 100644 --- a/bin/small_case.cc +++ b/bin/small_case.cc @@ -1,83 +1,83 @@ /** * @file small_case.cc * * @author Till Junge * * @date 12 Jan 2018 * * @brief small case for debugging * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/muSpectre_common.hh" #include "cell/cell_factory.hh" #include "materials/material_linear_elastic1.hh" #include "solver/solvers.hh" #include "solver/solver_cg.hh" #include #include using namespace muSpectre; int main() { constexpr Dim_t dim{twoD}; Ccoord_t resolution{11, 11}; Rcoord_t lengths{ - muGrid::CcoordOps::get_cube(11.)}; // {5.2e-9, 8.3e-9, 8.3e-9}; + muGrid::CcoordOps::get_cube(11.)}; // {5.2e-9, 8.3e-9, 8.3e-9}; Formulation form{Formulation::finite_strain}; auto rve{make_cell(resolution, lengths, form)}; auto & hard{MaterialLinearElastic1::make(rve, "hard", 210., .33)}; auto & soft{MaterialLinearElastic1::make(rve, "soft", 70., .33)}; for (auto && tup : akantu::enumerate(rve)) { auto & i = std::get<0>(tup); auto & pixel = std::get<1>(tup); if (i < 3) { hard.add_pixel(pixel); } else { soft.add_pixel(pixel); } } rve.initialise(); Real tol{1e-6}; Eigen::MatrixXd Del0{}; Del0 << 0, .1, 0, 0; Uint maxiter{31}; Dim_t verbose{3}; SolverCG cg{rve, tol, maxiter, static_cast(verbose)}; auto res = de_geus(rve, Del0, cg, tol, verbose); std::cout << res.grad.transpose() << std::endl; return 0; } diff --git a/language_bindings/python/bind_py_cell.cc b/language_bindings/python/bind_py_cell.cc index 78e1ff2..801e06c 100644 --- a/language_bindings/python/bind_py_cell.cc +++ b/language_bindings/python/bind_py_cell.cc @@ -1,262 +1,261 @@ /** * @file bind_py_cell.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief Python bindings for the cell factory function * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/muSpectre_common.hh" #include #include "cell/cell_factory.hh" #include "cell/cell_base.hh" #ifdef WITH_FFTWMPI #include "projection/fftwmpi_engine.hh" #endif #ifdef WITH_PFFT #include "projection/pfft_engine.hh" #endif #include #include #include #include "pybind11/eigen.h" #include #include using muSpectre::Ccoord_t; using muSpectre::Dim_t; using muSpectre::Formulation; using muSpectre::Rcoord_t; using pybind11::literals::operator""_a; namespace py = pybind11; /** * cell factory for specific FFT engine */ #ifdef WITH_MPI template void add_parallel_cell_factory_helper(py::module & mod, const char * name) { using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; mod.def(name, [](Ccoord res, Rcoord lens, Formulation form, size_t comm) { return make_parallel_cell, FFTEngine>( std::move(res), std::move(lens), std::move(form), std::move(Communicator(MPI_Comm(comm)))); }, "resolutions"_a, "lengths"_a = CcoordOps::get_cube(1.), "formulation"_a = Formulation::finite_strain, "communicator"_a = size_t(MPI_COMM_SELF)); } #endif /** * the cell factory is only bound for default template parameters */ template void add_cell_factory_helper(py::module & mod) { using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; mod.def("CellFactory", [](Ccoord res, Rcoord lens, Formulation form) { return make_cell(std::move(res), std::move(lens), std::move(form)); }, - "resolutions"_a, - "lengths"_a = muSpectre::CcoordOps::get_cube(1.), + "resolutions"_a, "lengths"_a = muGrid::CcoordOps::get_cube(1.), "formulation"_a = Formulation::finite_strain); #ifdef WITH_FFTWMPI - add_parallel_cell_factory_helper>( + add_parallel_cell_factory_helper>( mod, "FFTWMPICellFactory"); #endif #ifdef WITH_PFFT - add_parallel_cell_factory_helper>(mod, - "PFFTCellFactory"); + add_parallel_cell_factory_helper>( + mod, "PFFTCellFactory"); #endif } void add_cell_factory(py::module & mod) { - add_cell_factory_helper(mod); - add_cell_factory_helper(mod); + add_cell_factory_helper(mod); + add_cell_factory_helper(mod); } /** * CellBase for which the material and spatial dimension are identical */ template void add_cell_base_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "CellBase" << dim << 'd'; const std::string name = name_stream.str(); using sys_t = muSpectre::CellBase; py::class_(mod, name.c_str()) .def("__len__", &sys_t::size) .def("__iter__", [](sys_t & s) { return py::make_iterator(s.begin(), s.end()); }) .def("initialise", &sys_t::initialise, - "flags"_a = muSpectre::FFT_PlanFlags::estimate) + "flags"_a = muFFT::FFT_PlanFlags::estimate) .def( "directional_stiffness", [](sys_t & cell, py::EigenDRef & v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim * dim)) { std::stringstream err{}; err << "need array of shape (" << dim * dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } if (!cell.is_initialised()) { cell.initialise(); } const std::string out_name{"temp output for directional stiffness"}; const std::string in_name{"temp input for directional stiffness"}; constexpr bool create_tangent{true}; auto & K = cell.get_tangent(create_tangent); auto & input = cell.get_managed_T2_field(in_name); auto & output = cell.get_managed_T2_field(out_name); input.eigen() = v; cell.directional_stiffness(K, input, output); return output.eigen(); }, "δF"_a) .def("project", [](sys_t & cell, py::EigenDRef & v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim * dim)) { std::stringstream err{}; err << "need array of shape (" << dim * dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } if (!cell.is_initialised()) { cell.initialise(); } const std::string in_name{"temp input for projection"}; auto & input = cell.get_managed_T2_field(in_name); input.eigen() = v; cell.project(input); return input.eigen(); }, "field"_a) .def("get_strain", [](sys_t & s) { return s.get_strain().eigen(); }, py::return_value_policy::reference_internal) .def("get_stress", [](sys_t & s) { return Eigen::ArrayXXd(s.get_stress().eigen()); }) .def_property_readonly("size", &sys_t::size) .def("evaluate_stress_tangent", [](sys_t & cell, py::EigenDRef & v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim * dim)) { std::stringstream err{}; err << "need array of shape (" << dim * dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } auto & strain{cell.get_strain()}; strain.eigen() = v; auto stress_tgt{cell.evaluate_stress_tangent(strain)}; return std::tuple( std::get<0>(stress_tgt).eigen(), std::get<1>(stress_tgt).eigen()); }, "strain"_a) .def("evaluate_stress", [](sys_t & cell, py::EigenDRef & v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim * dim)) { std::stringstream err{}; err << "need array of shape (" << dim * dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } auto & strain{cell.get_strain()}; strain.eigen() = v; return cell.evaluate_stress(); }, "strain"_a, py::return_value_policy::reference_internal) .def("get_projection", &sys_t::get_projection) .def("get_subdomain_resolutions", &sys_t::get_subdomain_resolutions) .def("get_subdomain_locations", &sys_t::get_subdomain_locations) .def("get_domain_resolutions", &sys_t::get_domain_resolutions) .def("get_domain_lengths", &sys_t::get_domain_resolutions) .def("set_uniform_strain", [](sys_t & cell, py::EigenDRef & v) -> void { cell.set_uniform_strain(v); }, "strain"_a) .def("save_history_variables", &sys_t::save_history_variables); } void add_cell_base(py::module & mod) { py::class_(mod, "Cell") .def("get_globalised_internal_real_array", &muSpectre::Cell::get_globalised_internal_real_array, "unique_name"_a, "Convenience function to copy local (internal) fields of " "materials into a global field. At least one of the materials in " "the cell needs to contain an internal field named " "`unique_name`. If multiple materials contain such a field, they " "all need to be of same scalar type and same number of " "components. This does not work for split pixel cells or " "laminate pixel cells, as they can have multiple entries for the " "same pixel. Pixels for which no field named `unique_name` " "exists get an array of zeros." "\n" "Parameters:\n" "unique_name: fieldname to fill the global field with. At " "least one material must have such a field, or an " "Exception is raised.") .def("get_globalised_current_real_array", &muSpectre::Cell::get_globalised_current_real_array, "unique_name"_a) .def("get_globalised_old_real_array", &muSpectre::Cell::get_globalised_old_real_array, "unique_name"_a, "nb_steps_ago"_a = 1) .def("get_managed_real_array", &muSpectre::Cell::get_managed_real_array, "unique_name"_a, "nb_components"_a, "returns a field or nb_components real numbers per pixel"); add_cell_base_helper(mod); add_cell_base_helper(mod); } void add_cell(py::module & mod) { add_cell_factory(mod); auto cell{mod.def_submodule("cell")}; cell.doc() = "bindings for cells and cell factories"; add_cell_base(cell); } diff --git a/language_bindings/python/bind_py_common.cc b/language_bindings/python/bind_py_common.cc index fad517f..915fc27 100644 --- a/language_bindings/python/bind_py_common.cc +++ b/language_bindings/python/bind_py_common.cc @@ -1,174 +1,172 @@ /** * @file bind_py_common.cc * * @author Till Junge * * @date 08 Jan 2018 * * @brief Python bindings for the common part of µSpectre * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/muSpectre_common.hh" #include #include #include #include using muSpectre::Dim_t; using muSpectre::Real; using muSpectre::StressMeasure; using muSpectre::StrainMeasure; using muSpectre::Formulation; using pybind11::literals::operator""_a; namespace py = pybind11; template void add_get_cube_helper(py::module & mod) { std::stringstream name{}; name << "get_" << dim << "d_cube"; - mod.def(name.str().c_str(), &muSpectre::CcoordOps::get_cube, "size"_a, + mod.def(name.str().c_str(), &muGrid::CcoordOps::get_cube, "size"_a, "return a Ccoord with the value 'size' repeated in each dimension"); } template void add_get_hermitian_helper(py::module & mod) { - mod.def("get_hermitian_sizes", - &muSpectre::CcoordOps::get_hermitian_sizes, "full_sizes"_a, + mod.def("get_hermitian_sizes", &muGrid::CcoordOps::get_hermitian_sizes, + "full_sizes"_a, "return the hermitian sizes corresponding to the true sizes"); } template void add_get_ccoord_helper(py::module & mod) { - using Ccoord = muSpectre::Ccoord_t; + using Ccoord = muGrid::Ccoord_t; mod.def( "get_domain_ccoord", [](Ccoord resolutions, Dim_t index) { - return muSpectre::CcoordOps::get_ccoord(resolutions, Ccoord{}, - index); + return muGrid::CcoordOps::get_ccoord(resolutions, Ccoord{}, index); }, "resolutions"_a, "i"_a, "return the cell coordinate corresponding to the i'th cell in a grid of " "shape resolutions"); } void add_get_cube(py::module & mod) { add_get_cube_helper(mod); add_get_cube_helper(mod); add_get_cube_helper(mod); add_get_cube_helper(mod); add_get_hermitian_helper(mod); add_get_hermitian_helper(mod); add_get_ccoord_helper(mod); add_get_ccoord_helper(mod); } template void add_get_index_helper(py::module & mod) { - using Ccoord = muSpectre::Ccoord_t; + using Ccoord = muGrid::Ccoord_t; mod.def("get_domain_index", [](Ccoord sizes, Ccoord ccoord) { - return muSpectre::CcoordOps::get_index(sizes, Ccoord{}, - ccoord); + return muGrid::CcoordOps::get_index(sizes, Ccoord{}, ccoord); }, "sizes"_a, "ccoord"_a, "return the linear index corresponding to grid point 'ccoord' in a " "grid of size 'sizes'"); } void add_get_index(py::module & mod) { add_get_index_helper(mod); add_get_index_helper(mod); } template void add_Pixels_helper(py::module & mod) { std::stringstream name{}; name << "Pixels" << dim << "d"; - using Ccoord = muSpectre::Ccoord_t; - py::class_> Pixels(mod, name.str().c_str()); + using Ccoord = muGrid::Ccoord_t; + py::class_> Pixels(mod, name.str().c_str()); Pixels.def(py::init()); } void add_Pixels(py::module & mod) { - add_Pixels_helper(mod); - add_Pixels_helper(mod); + add_Pixels_helper(mod); + add_Pixels_helper(mod); } void add_common(py::module & mod) { py::enum_(mod, "Formulation") .value("finite_strain", Formulation::finite_strain) .value("small_strain", Formulation::small_strain); py::enum_(mod, "StressMeasure") .value("Cauchy", StressMeasure::Cauchy) .value("PK1", StressMeasure::PK1) .value("PK2", StressMeasure::PK2) .value("Kirchhoff", StressMeasure::Kirchhoff) .value("Biot", StressMeasure::Biot) .value("Mandel", StressMeasure::Mandel) .value("no_stress_", StressMeasure::no_stress_); py::enum_(mod, "StrainMeasure") .value("Gradient", StrainMeasure::Gradient) .value("Infinitesimal", StrainMeasure::Infinitesimal) .value("GreenLagrange", StrainMeasure::GreenLagrange) .value("Biot", StrainMeasure::Biot) .value("Log", StrainMeasure::Log) .value("Almansi", StrainMeasure::Almansi) .value("RCauchyGreen", StrainMeasure::RCauchyGreen) .value("LCauchyGreen", StrainMeasure::LCauchyGreen) .value("no_strain_", StrainMeasure::no_strain_); - py::enum_(mod, "FFT_PlanFlags") - .value("estimate", muSpectre::FFT_PlanFlags::estimate) - .value("measure", muSpectre::FFT_PlanFlags::measure) - .value("patient", muSpectre::FFT_PlanFlags::patient); + py::enum_(mod, "FFT_PlanFlags") + .value("estimate", muFFT::FFT_PlanFlags::estimate) + .value("measure", muFFT::FFT_PlanFlags::measure) + .value("patient", muFFT::FFT_PlanFlags::patient); py::enum_( mod, "FiniteDiff", "Distinguishes between different options of numerical differentiation;\n " " 1) 'forward' finite differences: ∂f/∂x ≈ (f(x+Δx) - f(x))/Δx\n 2) " "'backward' finite differences: ∂f/∂x ≈ (f(x) - f(x-Δx))/Δx\n 3) " "'centred' finite differences: ∂f/∂x ≈ (f(x+Δx) - f(x-Δx))/2Δx") .value("forward", muSpectre::FiniteDiff::forward) .value("backward", muSpectre::FiniteDiff::backward) .value("centred", muSpectre::FiniteDiff::centred); mod.def("banner", &muSpectre::banner, "name"_a, "year"_a, "copyright_holder"_a); add_get_cube(mod); add_Pixels(mod); add_get_index(mod); } diff --git a/language_bindings/python/bind_py_fftengine.cc b/language_bindings/python/bind_py_fftengine.cc index 6de828b..5145ec4 100644 --- a/language_bindings/python/bind_py_fftengine.cc +++ b/language_bindings/python/bind_py_fftengine.cc @@ -1,122 +1,127 @@ /** * @file bind_py_fftengine.cc * * @author Till Junge * * @date 17 Jan 2018 * * @brief Python bindings for the FFT engines * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ -#include "projection/fftw_engine.hh" +#include "bind_py_declarations.hh" + +#include #ifdef WITH_FFTWMPI -#include "projection/fftwmpi_engine.hh" +#include #endif #ifdef WITH_PFFT -#include "projection/pfft_engine.hh" +#include #endif -#include "bind_py_declarations.hh" #include #include #include -using muSpectre::Ccoord_t; -using muSpectre::Complex; -using muSpectre::Dim_t; +using muGrid::Ccoord_t; +using muGrid::Complex; +using muGrid::Dim_t; using pybind11::literals::operator""_a; namespace py = pybind11; template void add_engine_helper(py::module & mod, std::string name) { using Ccoord = Ccoord_t; using ArrayXXc = Eigen::Array; py::class_(mod, name.c_str()) #ifdef WITH_MPI .def(py::init([](Ccoord res, Dim_t nb_components, size_t comm) { return new Engine(res, nb_components, std::move(Communicator(MPI_Comm(comm)))); }), "resolutions"_a, "nb_components"_a, "communicator"_a = size_t(MPI_COMM_SELF)) #else .def(py::init()) #endif .def("fft", [](Engine & eng, py::EigenDRef v) { using Coll_t = typename Engine::GFieldCollection_t; using Field_t = typename Engine::Field_t; Coll_t coll{}; coll.initialise(eng.get_subdomain_resolutions(), eng.get_subdomain_locations()); - Field_t & temp{muSpectre::make_field( + Field_t & temp{muGrid::make_field( "temp_field", coll, eng.get_nb_components())}; temp.eigen() = v; return ArrayXXc{eng.fft(temp).eigen()}; }, "array"_a) .def("ifft", [](Engine & eng, py::EigenDRef v) { using Coll_t = typename Engine::GFieldCollection_t; using Field_t = typename Engine::Field_t; Coll_t coll{}; coll.initialise(eng.get_subdomain_resolutions(), eng.get_subdomain_locations()); - Field_t & temp{muSpectre::make_field( + Field_t & temp{muGrid::make_field( "temp_field", coll, eng.get_nb_components())}; eng.get_work_space().eigen() = v; eng.ifft(temp); return Eigen::ArrayXXd{temp.eigen()}; }, "array"_a) .def("initialise", &Engine::initialise, - "flags"_a = muSpectre::FFT_PlanFlags::estimate) + "flags"_a = muFFT::FFT_PlanFlags::estimate) .def("normalisation", &Engine::normalisation) .def("get_subdomain_resolutions", &Engine::get_subdomain_resolutions) .def("get_subdomain_locations", &Engine::get_subdomain_locations) .def("get_fourier_resolutions", &Engine::get_fourier_resolutions) .def("get_fourier_locations", &Engine::get_fourier_locations) .def("get_domain_resolutions", &Engine::get_domain_resolutions); } void add_fft_engines(py::module & mod) { auto fft{mod.def_submodule("fft")}; fft.doc() = "bindings for µSpectre's fft engines"; - add_engine_helper, muSpectre::twoD>( - fft, "FFTW_2d"); - add_engine_helper, - muSpectre::threeD>(fft, "FFTW_3d"); + add_engine_helper, muGrid::twoD>(fft, + "FFTW_2d"); + add_engine_helper, muGrid::threeD>( + fft, "FFTW_3d"); #ifdef WITH_FFTWMPI - add_engine_helper, twoD>(fft, "FFTWMPI_2d"); - add_engine_helper, threeD>(fft, "FFTWMPI_3d"); + add_engine_helper, muGrid::twoD>( + fft, "FFTWMPI_2d"); + add_engine_helper, muGrid::threeD>( + fft, "FFTWMPI_3d"); #endif #ifdef WITH_PFFT - add_engine_helper, twoD>(fft, "PFFT_2d"); - add_engine_helper, threeD>(fft, "PFFT_3d"); + add_engine_helper, muGrid::twoD>(fft, + "PFFT_2d"); + add_engine_helper, muGrid::threeD>( + fft, "PFFT_3d"); #endif add_projections(fft); } diff --git a/language_bindings/python/bind_py_field_collection.cc b/language_bindings/python/bind_py_field_collection.cc index 29924f7..3a5b81e 100644 --- a/language_bindings/python/bind_py_field_collection.cc +++ b/language_bindings/python/bind_py_field_collection.cc @@ -1,276 +1,276 @@ /** * file bind_py_field_collection.cc * * @author Till Junge * * @date 05 Jul 2018 * * @brief Python bindings for µSpectre field collections * * @section LICENSE * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/muSpectre_common.hh" #include #include #include #include #include #include #include using muSpectre::Complex; using muSpectre::Dim_t; using muSpectre::Int; using muSpectre::Real; using muSpectre::Uint; using pybind11::literals::operator""_a; namespace py = pybind11; template void add_field_collection(py::module & mod) { std::stringstream name_stream{}; name_stream << "_" << (FieldCollectionDerived::Global ? "Global" : "Local") << "FieldCollection_" << Dim << 'd'; const auto name{name_stream.str()}; - using FC_t = muSpectre::FieldCollectionBase; + using FC_t = muGrid::FieldCollectionBase; py::class_(mod, name.c_str()) .def("get_real_field", [](FC_t & field_collection, std::string unique_name) -> typename FC_t::template TypedField_t & { return field_collection.template get_typed_field( unique_name); }, "unique_name"_a, py::return_value_policy::reference_internal) .def("get_int_field", [](FC_t & field_collection, std::string unique_name) -> typename FC_t::template TypedField_t & { return field_collection.template get_typed_field(unique_name); }, "unique_name"_a, py::return_value_policy::reference_internal) .def("get_uint_field", [](FC_t & field_collection, std::string unique_name) -> typename FC_t::template TypedField_t & { return field_collection.template get_typed_field( unique_name); }, "unique_name"_a, py::return_value_policy::reference_internal) .def("get_complex_field", [](FC_t & field_collection, std::string unique_name) -> typename FC_t::template TypedField_t & { return field_collection.template get_typed_field( unique_name); }, "unique_name"_a, py::return_value_policy::reference_internal) .def("get_real_statefield", [](FC_t & field_collection, std::string unique_name) -> typename FC_t::template TypedStateField_t & { return field_collection.template get_typed_statefield( unique_name); }, "unique_name"_a, py::return_value_policy::reference_internal) .def("get_int_statefield", [](FC_t & field_collection, std::string unique_name) -> typename FC_t::template TypedStateField_t & { return field_collection.template get_typed_statefield( unique_name); }, "unique_name"_a, py::return_value_policy::reference_internal) .def("get_uint_statefield", [](FC_t & field_collection, std::string unique_name) -> typename FC_t::template TypedStateField_t & { return field_collection.template get_typed_statefield( unique_name); }, "unique_name"_a, py::return_value_policy::reference_internal) .def("get_complex_statefield", [](FC_t & field_collection, std::string unique_name) -> typename FC_t::template TypedStateField_t & { return field_collection.template get_typed_statefield( unique_name); }, "unique_name"_a, py::return_value_policy::reference_internal) .def_property_readonly( "field_names", &FC_t::get_field_names, "returns the names of all fields in this collection") .def_property_readonly("statefield_names", &FC_t::get_statefield_names, "returns the names of all state fields in this " "collection"); } namespace internal_fill { /** * Needed for static switch when adding fillers to fields (static * switch on whether they are global). Default case is for global * fields. */ template struct FillHelper { static void add_fill(PyField & py_field) { py_field.def( "fill_from_local", [](Field & field, const typename Field::LocalField_t & local) { field.fill_from_local(local); }, "local"_a, "Fills the content of a local field into a global field " "(modifies only the pixels that are not present in the " "local field"); } }; /** * Specialisation for local fields */ template struct FillHelper { static void add_fill(PyField & py_field) { py_field.def( "fill_from_global", [](Field & field, const typename Field::GlobalField_t & global) { field.fill_from_global(global); }, "global"_a, "Fills the content of a global field into a local field."); } }; } // namespace internal_fill template void add_field(py::module & mod, std::string dtype_name) { - using Field_t = muSpectre::TypedField; + using Field_t = muGrid::TypedField; std::stringstream name_stream{}; name_stream << (FieldCollection::Global ? "Global" : "Local") << "Field" << dtype_name << "_" << FieldCollection::spatial_dim(); std::string name{name_stream.str()}; using Ref_t = py::EigenDRef>; py::class_ py_field(mod, name.c_str()); py_field .def_property("array", [](Field_t & field) { return field.eigen(); }, [](Field_t & field, Ref_t mat) { field.eigen() = mat; }, "array of stored data") .def_property_readonly( "array", [](const Field_t & field) { return field.eigen(); }, "array of stored data") .def_property("vector", [](Field_t & field) { return field.eigenvec(); }, [](Field_t & field, Ref_t mat) { field.eigen() = mat; }, "flattened array of stored data") .def_property_readonly( "vector", [](const Field_t & field) { return field.eigenvec(); }, "flattened array of stored data"); using FillHelper_t = internal_fill::FillHelper; FillHelper_t::add_fill(py_field); } template void add_field_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << (FieldCollection::Global ? "Global" : "Local") << "Field" << "_" << Dim; std::string name{name_stream.str()}; - using Field_t = muSpectre::internal::FieldBase; + using Field_t = muGrid::internal::FieldBase; py::class_(mod, name.c_str()) .def_property_readonly("name", &Field_t::get_name, "field name") .def_property_readonly("collection", &Field_t::get_collection, "Collection containing this field") .def_property_readonly("nb_components", &Field_t::get_nb_components, "number of scalars stored per pixel in this field") .def_property_readonly("stored_type", [](const Field_t & field) { return field.get_stored_typeid().name(); }, "fundamental type of scalars stored in this field") .def_property_readonly("size", &Field_t::size, "number of pixels in this field") .def("set_zero", &Field_t::set_zero, "Set all components in the field to zero"); add_field(mod, "Real"); add_field(mod, "Int"); } template void add_statefield(py::module & mod, std::string dtype_name) { - using StateField_t = muSpectre::TypedStateField; + using StateField_t = muGrid::TypedStateField; std::stringstream name_stream{}; name_stream << (FieldCollection::Global ? "Global" : "Local") << "StateField" << dtype_name << "_" << FieldCollection::spatial_dim(); std::string name{name_stream.str()}; py::class_(mod, name.c_str()) .def_property_readonly("current_field", &StateField_t::get_current_field, "returns the current field value") .def("get_old_field", &StateField_t::get_old_field, "nb_steps_ago"_a = 1, "returns the value this field held 'nb_steps_ago' steps ago", py::return_value_policy::reference_internal); } template void add_statefield_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << (FieldCollection::Global ? "Global" : "Local") << "StateField" << "_" << Dim; std::string name{name_stream.str()}; - using StateField_t = muSpectre::StateFieldBase; + using StateField_t = muGrid::StateFieldBase; py::class_(mod, name.c_str()) .def_property_readonly("prefix", &StateField_t::get_prefix, "state field prefix") .def_property_readonly("collection", &StateField_t::get_collection, "Collection containing this field") .def_property_readonly("nb_memory", &StateField_t::get_nb_memory, "number of old states stored") .def_property_readonly( "stored_type", [](const StateField_t & field) { return field.get_stored_typeid().name(); }, "fundamental type of scalars stored in this field"); add_statefield(mod, "Real"); add_statefield(mod, "Int"); } template void add_field_collection_helper(py::module & mod) { - add_field_helper>(mod); - add_field_helper>(mod); + add_field_helper>(mod); + add_field_helper>(mod); - add_statefield_helper>(mod); - add_statefield_helper>(mod); + add_statefield_helper>(mod); + add_statefield_helper>(mod); - add_field_collection>(mod); - add_field_collection>(mod); + add_field_collection>(mod); + add_field_collection>(mod); } void add_field_collections(py::module & mod) { - add_field_collection_helper(mod); - add_field_collection_helper(mod); + add_field_collection_helper(mod); + add_field_collection_helper(mod); } diff --git a/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc index e7f7516..863b019 100644 --- a/language_bindings/python/bind_py_material.cc +++ b/language_bindings/python/bind_py_material.cc @@ -1,237 +1,235 @@ /** * @file bind_py_material.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief python bindings for µSpectre's materials * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/muSpectre_common.hh" #include "materials/material_base.hh" #include "materials/material_evaluator.hh" #include #include #include #include #include using muSpectre::Dim_t; using muSpectre::Real; using pybind11::literals::operator""_a; namespace py = pybind11; /* ---------------------------------------------------------------------- */ template void add_material_linear_elastic_generic1_helper(py::module & mod); /* ---------------------------------------------------------------------- */ template void add_material_linear_elastic_generic2_helper(py::module & mod); /** * python binding for the optionally objective form of Hooke's law */ template void add_material_linear_elastic1_helper(py::module & mod); template void add_material_linear_elastic2_helper(py::module & mod); template void add_material_linear_elastic3_helper(py::module & mod); template void add_material_linear_elastic4_helper(py::module & mod); template void add_material_hyper_elasto_plastic1_helper(py::module & mod); /* ---------------------------------------------------------------------- */ template class PyMaterialBase : public muSpectre::MaterialBase { public: /* Inherit the constructors */ using Parent = muSpectre::MaterialBase; using Parent::Parent; /* Trampoline (need one for each virtual function) */ void save_history_variables() override { PYBIND11_OVERLOAD_PURE(void, // Return type Parent, // Parent class save_history_variables); // Name of function in C++ // (must match Python name) } /* Trampoline (need one for each virtual function) */ void initialise() override { PYBIND11_OVERLOAD_PURE( void, // Return type Parent, // Parent class initialise); // Name of function in C++ (must match Python name) } void compute_stresses(const typename Parent::StrainField_t & F, typename Parent::StressField_t & P, muSpectre::Formulation form) override { PYBIND11_OVERLOAD_PURE( void, // Return type Parent, // Parent class compute_stresses, // Name of function in C++ (must match Python name) F, P, form); } void compute_stresses_tangent(const typename Parent::StrainField_t & F, typename Parent::StressField_t & P, typename Parent::TangentField_t & K, muSpectre::Formulation form) override { PYBIND11_OVERLOAD_PURE( void, /* Return type */ Parent, /* Parent class */ compute_stresses, /* Name of function in C++ (must match Python name) */ F, P, K, form); } }; template void add_material_evaluator(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialEvaluator_" << Dim << "d"; std::string name{name_stream.str()}; using MatEval_t = muSpectre::MaterialEvaluator; py::class_(mod, name.c_str()) .def(py::init>>()) .def("save_history_variables", &MatEval_t::save_history_variables, "for materials with state variables") .def("evaluate_stress", [](MatEval_t & mateval, py::EigenDRef & grad, muSpectre::Formulation form) { if ((grad.cols() != Dim) or (grad.rows() != Dim)) { std::stringstream err{}; err << "need matrix of shape (" << Dim << "×" << Dim << ") but got (" << grad.rows() << "×" << grad.cols() << ")."; throw std::runtime_error(err.str()); } return mateval.evaluate_stress(grad, form); }, "strain"_a, "formulation"_a, "Evaluates stress for a given strain and formulation " "(Piola-Kirchhoff 1 stress as a function of the placement gradient " "P = P(F) for formulation=Formulation.finite_strain and Cauchy " "stress as a function of the infinitesimal strain tensor σ = σ(ε) " "for formulation=Formulation.small_strain)") .def("evaluate_stress_tangent", [](MatEval_t & mateval, py::EigenDRef & grad, muSpectre::Formulation form) { if ((grad.cols() != Dim) or (grad.rows() != Dim)) { std::stringstream err{}; err << "need matrix of shape (" << Dim << "×" << Dim << ") but got (" << grad.rows() << "×" << grad.cols() << ")."; throw std::runtime_error(err.str()); } return mateval.evaluate_stress_tangent(grad, form); }, "strain"_a, "formulation"_a, "Evaluates stress and tangent moduli for a given strain and " "formulation (Piola-Kirchhoff 1 stress as a function of the " "placement gradient P = P(F) for " "formulation=Formulation.finite_strain and Cauchy stress as a " "function of the infinitesimal strain tensor σ = σ(ε) for " "formulation=Formulation.small_strain). The tangent moduli are K = " "∂P/∂F for formulation=Formulation.finite_strain and C = ∂σ/∂ε for " "formulation=Formulation.small_strain.") .def("estimate_tangent", [](MatEval_t & evaluator, py::EigenDRef & grad, muSpectre::Formulation form, const Real step, const muSpectre::FiniteDiff diff_type) { if ((grad.cols() != Dim) or (grad.rows() != Dim)) { std::stringstream err{}; err << "need matrix of shape (" << Dim << "×" << Dim << ") but got (" << grad.rows() << "×" << grad.cols() << ")."; throw std::runtime_error(err.str()); } return evaluator.estimate_tangent(grad, form, step, diff_type); }, "strain"_a, "formulation"_a, "Delta_x"_a, "difference_type"_a = muSpectre::FiniteDiff::centred, "Numerical estimate of the tangent modulus using finite " "differences. The finite difference scheme as well as the finite " "step size can be chosen. If there are no special circumstances, " "the default scheme of centred finite differences yields the most " "accurate results at an increased computational cost."); } template void add_material_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialBase_" << Dim << "d"; std::string name{name_stream.str()}; using Material = muSpectre::MaterialBase; using MaterialTrampoline = PyMaterialBase; - using FC_t = muSpectre::LocalFieldCollection; - using FCBase_t = muSpectre::FieldCollectionBase; + using FC_t = muGrid::LocalFieldCollection; + using FCBase_t = muGrid::FieldCollectionBase; py::class_>(mod, - name.c_str()) + std::shared_ptr>(mod, name.c_str()) .def(py::init()) .def("save_history_variables", &Material::save_history_variables) .def("list_fields", &Material::list_fields) .def("get_real_field", &Material::get_real_field, "field_name"_a, py::return_value_policy::reference_internal) .def("size", &Material::size) - .def("add_pixel", - [](Material & mat, muSpectre::Ccoord_t pix) { - mat.add_pixel(pix); - }, - "pixel"_a) + .def( + "add_pixel", + [](Material & mat, muGrid::Ccoord_t pix) { mat.add_pixel(pix); }, + "pixel"_a) .def_property_readonly("collection", [](Material & material) -> FCBase_t & { return material.get_collection(); }, "returns the field collection containing internal " "fields of this material"); add_material_linear_elastic1_helper(mod); add_material_linear_elastic2_helper(mod); add_material_linear_elastic3_helper(mod); add_material_linear_elastic4_helper(mod); add_material_hyper_elasto_plastic1_helper(mod); add_material_linear_elastic_generic1_helper(mod); add_material_linear_elastic_generic2_helper(mod); add_material_evaluator(mod); } void add_material(py::module & mod) { auto material{mod.def_submodule("material")}; material.doc() = "bindings for constitutive laws"; - add_material_helper(material); - add_material_helper(material); + add_material_helper(material); + add_material_helper(material); } diff --git a/language_bindings/python/bind_py_projections.cc b/language_bindings/python/bind_py_projections.cc index c10e55e..bb10282 100644 --- a/language_bindings/python/bind_py_projections.cc +++ b/language_bindings/python/bind_py_projections.cc @@ -1,191 +1,189 @@ /** * @file bind_py_projections.cc * * @author Till Junge * * @date 18 Jan 2018 * * @brief Python bindings for the Projection operators * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "projection/projection_small_strain.hh" #include "projection/projection_finite_strain.hh" #include "projection/projection_finite_strain_fast.hh" -#include "projection/fftw_engine.hh" +#include #ifdef WITH_FFTWMPI -#include "projection/fftwmpi_engine.hh" +#include #endif #ifdef WITH_PFFT -#include "projection/pfft_engine.hh" +#include #endif #include #include #include #include #include -using muSpectre::Dim_t; +using muGrid::Dim_t; using muSpectre::ProjectionBase; using pybind11::literals::operator""_a; namespace py = pybind11; /** * "Trampoline" class for handling the pure virtual methods, see * [http://pybind11.readthedocs.io/en/stable/advanced/classes.html#overriding-virtual-functions-in-python] * for details */ template class PyProjectionBase : public ProjectionBase { public: //! base class using Parent = ProjectionBase; //! field type on which projection is applied using Field_t = typename Parent::Field_t; void apply_projection(Field_t & field) override { PYBIND11_OVERLOAD_PURE(void, Parent, apply_projection, field); } Eigen::Map get_operator() override { PYBIND11_OVERLOAD_PURE(Eigen::Map, Parent, get_operator); } }; template void add_proj_helper(py::module & mod, std::string name_start) { - using Ccoord = muSpectre::Ccoord_t; - using Rcoord = muSpectre::Rcoord_t; + using Ccoord = muGrid::Ccoord_t; + using Rcoord = muGrid::Rcoord_t; using Field_t = typename Proj::Field_t; static_assert(DimS == DimM, "currently only for DimS==DimM"); std::stringstream name{}; name << name_start << '_' << DimS << 'd'; py::class_(mod, name.str().c_str()) #ifdef WITH_MPI .def(py::init([](Ccoord res, Rcoord lengths, const std::string & fft, size_t comm) { if (fft == "fftw") { - auto engine = std::make_unique>( + auto engine = std::make_unique>( res, Proj::NbComponents(), std::move(Communicator(MPI_Comm(comm)))); return Proj(std::move(engine), lengths); #else .def(py::init([](Ccoord res, Rcoord lengths, const std::string & fft) { if (fft == "fftw") { - auto engine = std::make_unique>( + auto engine = std::make_unique>( res, Proj::NbComponents()); return Proj(std::move(engine), lengths); #endif #ifdef WITH_FFTWMPI } else if (fft == "fftwmpi") { - auto engine = std::make_unique>( + auto engine = std::make_unique>( res, Proj::NbComponents(), std::move(Communicator(MPI_Comm(comm)))); return Proj(std::move(engine), lengths); #endif #ifdef WITH_PFFT } else if (fft == "pfft") { - auto engine = std::make_unique>( + auto engine = std::make_unique>( res, Proj::NbComponents(), std::move(Communicator(MPI_Comm(comm)))); return Proj(std::move(engine), lengths); #endif } else { throw std::runtime_error("Unknown FFT engine '" + fft + "' specified."); } }), "resolutions"_a, "lengths"_a, #ifdef WITH_MPI "fft"_a = "fftw", "communicator"_a = size_t(MPI_COMM_SELF)) #else "fft"_a = "fftw") #endif .def("initialise", &Proj::initialise, - "flags"_a = muSpectre::FFT_PlanFlags::estimate, + "flags"_a = muFFT::FFT_PlanFlags::estimate, "initialises the fft engine (plan the transform)") .def("apply_projection", [](Proj & proj, py::EigenDRef v) { - typename muSpectre::FFTEngineBase::GFieldCollection_t coll{}; - Eigen::Index subdomain_size = muSpectre::CcoordOps::get_size( - proj.get_subdomain_resolutions()); + typename muFFT::FFTEngineBase::GFieldCollection_t coll{}; + Eigen::Index subdomain_size = + muGrid::CcoordOps::get_size(proj.get_subdomain_resolutions()); if (v.rows() != DimS * DimM || v.cols() != subdomain_size) { throw std::runtime_error("Expected input array of shape (" + std::to_string(DimS * DimM) + ", " + std::to_string(subdomain_size) + "), but input array has shape (" + std::to_string(v.rows()) + ", " + std::to_string(v.cols()) + ")."); } coll.initialise(proj.get_subdomain_resolutions(), proj.get_subdomain_locations()); - Field_t & temp{muSpectre::make_field( + Field_t & temp{muGrid::make_field( "temp_field", coll, proj.get_nb_components())}; temp.eigen() = v; proj.apply_projection(temp); return Eigen::ArrayXXd{temp.eigen()}; }) .def("get_operator", &Proj::get_operator) .def( "get_formulation", &Proj::get_formulation, "return a Formulation enum indicating whether the projection is small" " or finite strain") .def("get_subdomain_resolutions", &Proj::get_subdomain_resolutions) .def("get_subdomain_locations", &Proj::get_subdomain_locations) .def("get_domain_resolutions", &Proj::get_domain_resolutions) .def("get_domain_lengths", &Proj::get_domain_resolutions); } void add_proj_dispatcher(py::module & mod) { + add_proj_helper, + muGrid::twoD>(mod, "ProjectionSmallStrain"); add_proj_helper< - muSpectre::ProjectionSmallStrain, - muSpectre::twoD>(mod, "ProjectionSmallStrain"); - add_proj_helper< - muSpectre::ProjectionSmallStrain, - muSpectre::threeD>(mod, "ProjectionSmallStrain"); + muSpectre::ProjectionSmallStrain, + muGrid::threeD>(mod, "ProjectionSmallStrain"); + add_proj_helper, + muGrid::twoD>(mod, "ProjectionFiniteStrain"); add_proj_helper< - muSpectre::ProjectionFiniteStrain, - muSpectre::twoD>(mod, "ProjectionFiniteStrain"); - add_proj_helper< - muSpectre::ProjectionFiniteStrain, - muSpectre::threeD>(mod, "ProjectionFiniteStrain"); + muSpectre::ProjectionFiniteStrain, + muGrid::threeD>(mod, "ProjectionFiniteStrain"); add_proj_helper< - muSpectre::ProjectionFiniteStrainFast, - muSpectre::twoD>(mod, "ProjectionFiniteStrainFast"); - add_proj_helper, - muSpectre::threeD>(mod, "ProjectionFiniteStrainFast"); + muSpectre::ProjectionFiniteStrainFast, + muGrid::twoD>(mod, "ProjectionFiniteStrainFast"); + add_proj_helper< + muSpectre::ProjectionFiniteStrainFast, + muGrid::threeD>(mod, "ProjectionFiniteStrainFast"); } void add_projections(py::module & mod) { add_proj_dispatcher(mod); } diff --git a/src/common/geometry.hh b/src/common/geometry.hh index cdee272..05857c4 100644 --- a/src/common/geometry.hh +++ b/src/common/geometry.hh @@ -1,263 +1,265 @@ /** * @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/muSpectre_common.hh" -#include "common/tensor_algebra.hh" +#include #include #include #include #include #ifndef SRC_COMMON_GEOMETRY_HH_ #define SRC_COMMON_GEOMETRY_HH_ 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}; }; } // namespace internal template ::value> class Rotator { public: 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"); using Angles_t = Eigen::Matrix; using RotMat_t = Eigen::Matrix; //! Default constructor Rotator() = delete; explicit Rotator(const Eigen::Ref & angles) : angles{angles}, rot_mat{this->compute_rotation_matrix()} {} //! Copy constructor Rotator(const Rotator & other) = default; //! Move constructor Rotator(Rotator && other) = default; //! Destructor virtual ~Rotator() = default; //! Copy assignment operator Rotator & operator=(const Rotator & other) = default; //! Move assignment operator Rotator & operator=(Rotator && other) = default; /** * Applies the rotation into the frame defined by 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 defined by 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 rot_mat; } protected: inline RotMat_t compute_rotation_matrix(); EIGEN_MAKE_ALIGNED_OPERATOR_NEW; Angles_t angles; RotMat_t rot_mat; private: }; namespace internal { template struct RotationMatrixComputer {}; template struct RotationMatrixComputer { constexpr static Dim_t Dim{twoD}; using RotMat_t = typename Rotator::RotMat_t; using Angles_t = typename Rotator::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 RotationMatrixComputer { constexpr static Dim_t Dim{threeD}; using RotMat_t = typename Rotator::RotMat_t; using Angles_t = typename Rotator::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."); } } } }; } // namespace internal /* ---------------------------------------------------------------------- */ template auto Rotator::compute_rotation_matrix() -> RotMat_t { return internal::RotationMatrixComputer::compute(this->angles); } 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}; + constexpr Dim_t Dim{muGrid::EigenCheck::tensor_dim::value}; auto && rotator_forward{ Matrices::outer_under(R.transpose(), R.transpose())}; auto && rotator_back = Matrices::outer_under(R, R); // Clarification. When I return this value as an // expression, clang segfaults or returns an uninitialised // tensor, hence the explicit cast into a T4Mat. - return T4Mat(rotator_back * input * rotator_forward); + return muGrid::T4Mat(rotator_back * input * rotator_forward); } }; } // namespace internal /* ---------------------------------------------------------------------- */ template template auto Rotator::rotate(In_t && input) -> decltype(auto) { - constexpr Dim_t tensor_rank{EigenCheck::tensor_rank::value}; + constexpr Dim_t tensor_rank{ + muGrid::EigenCheck::tensor_rank::value}; return internal::RotationHelper::rotate( std::forward(input), this->rot_mat); } /* ---------------------------------------------------------------------- */ template template auto Rotator::rotate_back(In_t && input) -> decltype(auto) { - constexpr Dim_t tensor_rank{EigenCheck::tensor_rank::value}; + constexpr Dim_t tensor_rank{ + muGrid::EigenCheck::tensor_rank::value}; return internal::RotationHelper::rotate( std::forward(input), this->rot_mat.transpose()); } } // namespace muSpectre #endif // SRC_COMMON_GEOMETRY_HH_ diff --git a/src/libmufft/mufft_common.hh b/src/libmufft/mufft_common.hh index f42fc1d..86266a6 100644 --- a/src/libmufft/mufft_common.hh +++ b/src/libmufft/mufft_common.hh @@ -1,66 +1,70 @@ /** * @file mufft_common.hh * * @author Till Junge * * @date 24 Jan 2019 * * @brief Small definitions of commonly used types throughout µFFT * * Copyright © 2019 Till Junge * * µFFT is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µFFT 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 Lesser General Public License * along with µFFT; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #ifndef SRC_LIBFFT_MUFFT_COMMON_HH_ #define SRC_LIBFFT_MUFFT_COMMON_HH_ namespace muFFT { using muGrid::Dim_t; using muGrid::Complex; using muGrid::Int; using muGrid::Real; using muGrid::Uint; using muGrid::Ccoord_t; using muGrid::Rcoord_t; using muGrid::optional; + using muGrid::oneD; + using muGrid::threeD; + using muGrid::twoD; + /** * Planner flags for FFT (follows FFTW, hopefully this choice will * be compatible with alternative FFT implementations) * @enum muFFT::FFT_PlanFlags */ enum class FFT_PlanFlags { estimate, //!< cheapest plan for slowest execution measure, //!< more expensive plan for fast execution patient //!< very expensive plan for fastest execution }; } // namespace muFFT #endif // SRC_LIBFFT_MUFFT_COMMON_HH_ diff --git a/src/materials/material_base.cc b/src/materials/material_base.cc index dc6ac01..6b66e38 100644 --- a/src/materials/material_base.cc +++ b/src/materials/material_base.cc @@ -1,108 +1,109 @@ /** * @file material_base.cc * * @author Till Junge * * @date 01 Nov 2017 * * @brief implementation of materi * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "materials/material_base.hh" namespace muSpectre { //----------------------------------------------------------------------------// template MaterialBase::MaterialBase(std::string name) : name(name) { static_assert((DimM == oneD) || (DimM == twoD) || (DimM == threeD), "only 1, 2, or threeD supported"); } /* ---------------------------------------------------------------------- */ template const std::string & MaterialBase::get_name() const { return this->name; } /* ---------------------------------------------------------------------- */ template void MaterialBase::add_pixel(const Ccoord & ccoord) { this->internal_fields.add_pixel(ccoord); } /* ---------------------------------------------------------------------- */ template void MaterialBase::compute_stresses(const Field_t & F, Field_t & P, Formulation form) { this->compute_stresses(StrainField_t::check_ref(F), StressField_t::check_ref(P), form); } /* ---------------------------------------------------------------------- */ template void MaterialBase::compute_stresses_tangent(const Field_t & F, Field_t & P, Field_t & K, Formulation form) { this->compute_stresses_tangent(StrainField_t::check_ref(F), StressField_t::check_ref(P), TangentField_t::check_ref(K), form); } //----------------------------------------------------------------------------// template auto MaterialBase::get_real_field(std::string field_name) -> EigenMap { if (not this->internal_fields.check_field_exists(field_name)) { std::stringstream err{}; err << "Field '" << field_name << "' does not exist in material '" << this->name << "'."; throw muGrid::FieldCollectionError(err.str()); } auto & field{this->internal_fields[field_name]}; if (field.get_stored_typeid().hash_code() != typeid(Real).hash_code()) { std::stringstream err{}; err << "Field '" << field_name << "' is not real-valued"; throw muGrid::FieldCollectionError(err.str()); } - return static_cast &>(field).eigen(); + return static_cast &>(field) + .eigen(); } /* ---------------------------------------------------------------------- */ template std::vector MaterialBase::list_fields() const { return this->internal_fields.list_fields(); } template class MaterialBase<2, 2>; template class MaterialBase<2, 3>; template class MaterialBase<3, 3>; } // namespace muSpectre diff --git a/src/materials/material_linear_elastic_generic2.cc b/src/materials/material_linear_elastic_generic2.cc index 5dfa4c6..cce7a7f 100644 --- a/src/materials/material_linear_elastic_generic2.cc +++ b/src/materials/material_linear_elastic_generic2.cc @@ -1,68 +1,69 @@ /** * @file material_linear_elastic_generic2.cc * * @author Till Junge * * @date 20 Dec 2018 * * @brief Implementation for generic linear elastic law with eigenstrains * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "material_linear_elastic_generic2.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ template MaterialLinearElasticGeneric2::MaterialLinearElasticGeneric2( const std::string & name, const CInput_t & C_voigt) : Parent{name}, worker{name, C_voigt}, - eigen_field{muGrid::make_field("Eigenstrain", this->internal_fields)}, + eigen_field{ + muGrid::make_field("Eigenstrain", this->internal_fields)}, internal_variables(eigen_field.get_const_map()) {} /* ---------------------------------------------------------------------- */ template void MaterialLinearElasticGeneric2::add_pixel( const Ccoord_t & /*pixel*/) { throw std::runtime_error("this material needs pixels with and eigenstrain"); } /* ---------------------------------------------------------------------- */ template void MaterialLinearElasticGeneric2::add_pixel( const Ccoord_t & pixel, const StrainTensor & E_eig) { this->internal_fields.add_pixel(pixel); Eigen::Map> strain_array( E_eig.data()); this->eigen_field.push_back(strain_array); } template class MaterialLinearElasticGeneric2; template class MaterialLinearElasticGeneric2; template class MaterialLinearElasticGeneric2; } // namespace muSpectre diff --git a/src/projection/projection_finite_strain.cc b/src/projection/projection_finite_strain.cc index 08faa4b..9a87192 100644 --- a/src/projection/projection_finite_strain.cc +++ b/src/projection/projection_finite_strain.cc @@ -1,97 +1,98 @@ /** * @file projection_finite_strain.cc * * @author Till Junge * * @date 05 Dec 2017 * * @brief implementation of standard finite strain projection operator * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "projection/projection_finite_strain.hh" #include #include #include #include #include "Eigen/Dense" namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionFiniteStrain::ProjectionFiniteStrain( FFTEngine_ptr engine, Rcoord lengths) : Parent{std::move(engine), lengths, Formulation::finite_strain} { for (auto res : this->fft_engine->get_domain_resolutions()) { if (res % 2 == 0) { throw ProjectionError( "Only an odd number of gridpoints in each direction is supported"); } } } /* ---------------------------------------------------------------------- */ template - void ProjectionFiniteStrain::initialise(muFFT::FFT_PlanFlags flags) { + void + ProjectionFiniteStrain::initialise(muFFT::FFT_PlanFlags flags) { Parent::initialise(flags); muFFT::FFT_freqs fft_freqs(this->fft_engine->get_domain_resolutions(), - this->domain_lengths); + this->domain_lengths); for (auto && tup : akantu::zip(*this->fft_engine, this->Ghat)) { const auto & ccoord = std::get<0>(tup); auto & G = std::get<1>(tup); auto xi = fft_freqs.get_unit_xi(ccoord); //! this is simplifiable using Curnier's Méthodes numériques, 6.69(c) G = Matrices::outer_under(Matrices::I2(), xi * xi.transpose()); // The commented block below corresponds to the original // definition of the operator in de Geus et // al. (https://doi.org/10.1016/j.cma.2016.12.032). However, // they use a bizarre definition of the double contraction // between fourth-order and second-order tensors that has a // built-in transpose operation (i.e., C = A:B <-> AᵢⱼₖₗBₗₖ = // Cᵢⱼ , note the inverted ₗₖ instead of ₖₗ), here, we define // the double contraction without the transposition. As a // result, the Projection operator produces the transpose of de // Geus's // for (Dim_t im = 0; im < DimS; ++im) { // for (Dim_t j = 0; j < DimS; ++j) { // for (Dim_t l = 0; l < DimS; ++l) { // get(G, im, j, l, im) = xi(j)*xi(l); // } // } // } } if (this->get_subdomain_locations() == Ccoord{}) { this->Ghat[0].setZero(); } } template class ProjectionFiniteStrain; template class ProjectionFiniteStrain; } // namespace muSpectre diff --git a/src/projection/projection_finite_strain_fast.cc b/src/projection/projection_finite_strain_fast.cc index 1766e08..b63350d 100644 --- a/src/projection/projection_finite_strain_fast.cc +++ b/src/projection/projection_finite_strain_fast.cc @@ -1,104 +1,105 @@ /** * @file projection_finite_strain_fast.cc * * @author Till Junge * * @date 12 Dec 2017 * * @brief implementation for fast projection in finite strain * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "projection/projection_finite_strain_fast.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionFiniteStrainFast::ProjectionFiniteStrainFast( FFTEngine_ptr engine, Rcoord lengths) : Parent{std::move(engine), lengths, Formulation::finite_strain}, xiField{muGrid::make_field("Projection Operator", - this->projection_container)}, + this->projection_container)}, xis(xiField) { for (auto res : this->fft_engine->get_domain_resolutions()) { if (res % 2 == 0) { throw ProjectionError( "Only an odd number of gridpoints in each direction is supported"); } } } /* ---------------------------------------------------------------------- */ template - void ProjectionFiniteStrainFast::initialise(muFFT::FFT_PlanFlags flags) { + void ProjectionFiniteStrainFast::initialise( + muFFT::FFT_PlanFlags flags) { Parent::initialise(flags); muFFT::FFT_freqs fft_freqs(this->fft_engine->get_domain_resolutions(), - this->domain_lengths); + this->domain_lengths); for (auto && tup : akantu::zip(*this->fft_engine, this->xis)) { const auto & ccoord = std::get<0>(tup); auto & xi = std::get<1>(tup); xi = fft_freqs.get_unit_xi(ccoord); } if (this->get_subdomain_locations() == Ccoord{}) { this->xis[0].setZero(); } } /* ---------------------------------------------------------------------- */ template void ProjectionFiniteStrainFast::apply_projection(Field_t & field) { Grad_map field_map{this->fft_engine->fft(field)}; Real factor = this->fft_engine->normalisation(); for (auto && tup : akantu::zip(this->xis, field_map)) { auto & xi{std::get<0>(tup)}; auto & f{std::get<1>(tup)}; f = factor * ((f * xi).eval() * xi.transpose()); } this->fft_engine->ifft(field); } /* ---------------------------------------------------------------------- */ template Eigen::Map ProjectionFiniteStrainFast::get_operator() { return this->xiField.dyn_eigen(); } /* ---------------------------------------------------------------------- */ template std::array ProjectionFiniteStrainFast::get_strain_shape() const { return std::array{DimM, DimM}; } /* ---------------------------------------------------------------------- */ template class ProjectionFiniteStrainFast; template class ProjectionFiniteStrainFast; } // namespace muSpectre diff --git a/src/projection/projection_small_strain.cc b/src/projection/projection_small_strain.cc index 5705537..89b20f7 100644 --- a/src/projection/projection_small_strain.cc +++ b/src/projection/projection_small_strain.cc @@ -1,90 +1,91 @@ /** * @file projection_small_strain.cc * * @author Till Junge * * @date 14 Jan 2018 * * @brief Implementation for ProjectionSmallStrain * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "projection/projection_small_strain.hh" #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionSmallStrain::ProjectionSmallStrain(FFTEngine_ptr engine, Rcoord lengths) : Parent{std::move(engine), lengths, Formulation::small_strain} { for (auto res : this->fft_engine->get_domain_resolutions()) { if (res % 2 == 0) { throw ProjectionError( "Only an odd number of gridpoints in each direction is supported"); } } } /* ---------------------------------------------------------------------- */ template - void ProjectionSmallStrain::initialise(muFFT::FFT_PlanFlags flags) { + void + ProjectionSmallStrain::initialise(muFFT::FFT_PlanFlags flags) { using muGrid::get; Parent::initialise(flags); - + muFFT::FFT_freqs fft_freqs(this->fft_engine->get_domain_resolutions(), - this->domain_lengths); + this->domain_lengths); for (auto && tup : akantu::zip(*this->fft_engine, this->Ghat)) { const auto & ccoord = std::get<0>(tup); auto & G = std::get<1>(tup); auto xi = fft_freqs.get_unit_xi(ccoord); auto kron = [](const Dim_t i, const Dim_t j) -> Real { return (i == j) ? 1. : 0.; }; for (Dim_t i{0}; i < DimS; ++i) { for (Dim_t j{0}; j < DimS; ++j) { for (Dim_t l{0}; l < DimS; ++l) { for (Dim_t m{0}; m < DimS; ++m) { Real & g = get(G, i, j, l, m); g = 0.5 * (xi(i) * kron(j, l) * xi(m) + xi(i) * kron(j, m) * xi(l) + xi(j) * kron(i, l) * xi(m) + xi(j) * kron(i, m) * xi(l)) - xi(i) * xi(j) * xi(l) * xi(m); } } } } } if (this->get_subdomain_locations() == Ccoord{}) { this->Ghat[0].setZero(); } } template class ProjectionSmallStrain; template class ProjectionSmallStrain; } // namespace muSpectre diff --git a/src/solver/deprecated_solver_cg.cc b/src/solver/deprecated_solver_cg.cc index a6bfb5f..dc92c93 100644 --- a/src/solver/deprecated_solver_cg.cc +++ b/src/solver/deprecated_solver_cg.cc @@ -1,139 +1,140 @@ /** * @file deprecated_solver_cg.cc * * @author Till Junge * * @date 20 Dec 2017 * * @brief Implementation of cg solver * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "solver/deprecated_solver_cg.hh" #include "solver/solver_common.hh" #include #include #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template DeprecatedSolverCG::DeprecatedSolverCG(Cell_t & cell, Real tol, Uint maxiter, bool verbose) : Parent(cell, tol, maxiter, verbose), r_k{muGrid::make_field( "residual r_k", this->collection)}, - p_k{muGrid::make_field("search direction r_k", this->collection)}, - Ap_k{muGrid::make_field("Effect of tangent A*p_k", this->collection)} { - } + p_k{muGrid::make_field("search direction r_k", + this->collection)}, + Ap_k{muGrid::make_field("Effect of tangent A*p_k", + this->collection)} {} /* ---------------------------------------------------------------------- */ template void DeprecatedSolverCG::solve(const Field_t & rhs, Field_t & x_f) { x_f.eigenvec() = this->solve(rhs.eigenvec(), x_f.eigenvec()); } //----------------------------------------------------------------------------// template typename DeprecatedSolverCG::SolvVectorOut DeprecatedSolverCG::solve(const SolvVectorInC rhs, SolvVectorIn x_0) { const muFFT::Communicator & comm = this->cell.get_communicator(); // Following implementation of algorithm 5.2 in Nocedal's Numerical // Optimization (p. 112) auto r = this->r_k.eigen(); auto p = this->p_k.eigen(); auto Ap = this->Ap_k.eigen(); typename Field_t::EigenMap_t x(x_0.data(), r.rows(), r.cols()); // initialisation of algo r = this->cell.directional_stiffness_with_copy(x); r -= typename Field_t::ConstEigenMap_t(rhs.data(), r.rows(), r.cols()); p = -r; this->converged = false; Real rdr = comm.sum((r * r).sum()); Real rhs_norm2 = comm.sum(rhs.squaredNorm()); Real tol2 = muGrid::ipow(this->tol, 2) * rhs_norm2; size_t count_width{}; // for output formatting in verbose case if (this->verbose) { count_width = size_t(std::log10(this->maxiter)) + 1; } for (Uint i = 0; i < this->maxiter && (rdr > tol2 || i == 0); ++i, ++this->counter) { Ap = this->cell.directional_stiffness_with_copy(p); Real alpha = rdr / comm.sum((p * Ap).sum()); x += alpha * p; r += alpha * Ap; Real new_rdr = comm.sum((r * r).sum()); Real beta = new_rdr / rdr; rdr = new_rdr; if (this->verbose && comm.rank() == 0) { std::cout << " at CG step " << std::setw(count_width) << i << ": |r|/|b| = " << std::setw(15) << std::sqrt(rdr / rhs_norm2) << ", cg_tol = " << this->tol << std::endl; } p = -r + beta * p; } if (rdr < tol2) { this->converged = true; } else { std::stringstream err{}; err << " After " << this->counter << " steps, the solver " << " FAILED with |r|/|b| = " << std::setw(15) << std::sqrt(rdr / rhs_norm2) << ", cg_tol = " << this->tol << std::endl; throw ConvergenceError("Conjugate gradient has not converged." + err.str()); } return x_0; } /* ---------------------------------------------------------------------- */ template typename DeprecatedSolverCG::Tg_req_t DeprecatedSolverCG::get_tangent_req() const { return tangent_requirement; } template class DeprecatedSolverCG; // template class DeprecatedSolverCG; template class DeprecatedSolverCG; } // namespace muSpectre diff --git a/tests/header_test_ccoord_operations.cc b/tests/header_test_ccoord_operations.cc index 9073e5f..65d6e3b 100644 --- a/tests/header_test_ccoord_operations.cc +++ b/tests/header_test_ccoord_operations.cc @@ -1,159 +1,156 @@ /** * @file test_ccoord_operations.cc * * @author Till Junge * * @date 03 Dec 2017 * * @brief tests for cell coordinate operations * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include #include "tests/test_goodies.hh" #include "tests.hh" -using namespace muSpectre; - namespace muGrid { - + using muSpectre::tol; BOOST_AUTO_TEST_SUITE(ccoords_operations); - BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_cube, Fix, - testGoodies::dimlist, Fix) { + BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_cube, Fix, testGoodies::dimlist, Fix) { constexpr auto dim{Fix::dim}; using Ccoord = Ccoord_t; constexpr Dim_t size{5}; constexpr Ccoord cube = CcoordOps::get_cube(size); Ccoord ref_cube; for (Dim_t i = 0; i < dim; ++i) { ref_cube[i] = size; } BOOST_CHECK_EQUAL_COLLECTIONS(ref_cube.begin(), ref_cube.end(), cube.begin(), cube.end()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_hermitian, Fix, testGoodies::dimlist, Fix) { constexpr auto dim{Fix::dim}; using Ccoord = Ccoord_t; constexpr Dim_t size{5}; constexpr Ccoord cube = CcoordOps::get_cube(size); constexpr Ccoord herm = CcoordOps::get_hermitian_sizes(cube); Ccoord ref_cube = cube; ref_cube.back() = (cube.back() + 1) / 2; BOOST_CHECK_EQUAL_COLLECTIONS(ref_cube.begin(), ref_cube.end(), herm.begin(), herm.end()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_size, Fix, testGoodies::dimlist, Fix) { constexpr auto dim{Fix::dim}; using Ccoord = Ccoord_t; constexpr Dim_t size{5}; constexpr Ccoord cube = CcoordOps::get_cube(size); BOOST_CHECK_EQUAL(CcoordOps::get_size(cube), ipow(size, dim)); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_stride_size, Fix, testGoodies::dimlist, Fix) { constexpr auto dim{Fix::dim}; using Ccoord = Ccoord_t; constexpr Dim_t size{5}; constexpr Ccoord cube = CcoordOps::get_cube(size); constexpr Ccoord stride = CcoordOps::get_default_strides(cube); BOOST_CHECK_EQUAL(CcoordOps::get_size_from_strides(cube, stride), ipow(size, dim)); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_index, Fix, testGoodies::dimlist, Fix) { constexpr auto dim{Fix::dim}; using Ccoord = Ccoord_t; testGoodies::RandRange rng; Ccoord sizes{}; for (Dim_t i{0}; i < dim; ++i) { sizes[i] = rng.randval(2, 5); } Ccoord stride = CcoordOps::get_default_strides(sizes); Ccoord locations{}; const size_t nb_pix{CcoordOps::get_size(sizes)}; for (size_t i{0}; i < nb_pix; ++i) { BOOST_CHECK_EQUAL( i, CcoordOps::get_index_from_strides( stride, CcoordOps::get_ccoord(sizes, locations, i))); } } BOOST_AUTO_TEST_CASE(vector_test) { constexpr Ccoord_t c3{1, 2, 3}; constexpr Ccoord_t c2{c3[0], c3[1]}; constexpr Rcoord_t s3{1.3, 2.8, 5.7}; constexpr Rcoord_t s2{s3[0], s3[1]}; Eigen::Matrix v2; v2 << s3[0], s3[1]; Eigen::Matrix v3; v3 << s3[0], s3[1], s3[2]; auto vec2{CcoordOps::get_vector(c2, v2(1))}; auto vec3{CcoordOps::get_vector(c3, v3(1))}; for (Dim_t i = 0; i < twoD; ++i) { BOOST_CHECK_EQUAL(c2[i] * v2(1), vec2[i]); } for (Dim_t i = 0; i < threeD; ++i) { BOOST_CHECK_EQUAL(c3[i] * v3(1), vec3[i]); } vec2 = CcoordOps::get_vector(c2, v2); vec3 = CcoordOps::get_vector(c3, v3); for (Dim_t i = 0; i < twoD; ++i) { BOOST_CHECK_EQUAL(c2[i] * v2(i), vec2[i]); BOOST_CHECK_EQUAL(vec2[i], CcoordOps::get_vector(c2, s2)[i]); } for (Dim_t i = 0; i < threeD; ++i) { BOOST_CHECK_EQUAL(c3[i] * v3(i), vec3[i]); BOOST_CHECK_EQUAL(vec3[i], CcoordOps::get_vector(c3, s3)[i]); } } BOOST_AUTO_TEST_SUITE_END(); } // namespace muGrid diff --git a/tests/header_test_eigen_tools.cc b/tests/header_test_eigen_tools.cc index 55270d2..06d638c 100644 --- a/tests/header_test_eigen_tools.cc +++ b/tests/header_test_eigen_tools.cc @@ -1,132 +1,133 @@ /** * @file header_test_eigen_tools.cc * * @author Till Junge * * @date 07 Mar 2018 * * @brief test the eigen_tools * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ -#include "common/eigen_tools.hh" #include "tests.hh" +#include +#include -namespace muSpectre { +namespace muGrid { BOOST_AUTO_TEST_SUITE(eigen_tools); BOOST_AUTO_TEST_CASE(exponential_test) { using Mat_t = Eigen::Matrix; Mat_t input{}; input << 0, .25 * pi, 0, .25 * pi, 0, 0, 0, 0, 1; Mat_t output{}; output << 1.32460909, 0.86867096, 0, 0.86867096, 1.32460909, 0, 0, 0, 2.71828183; auto my_output{expm(input)}; Real error{(my_output - output).norm()}; BOOST_CHECK_LT(error, 1e-8); if (error >= 1e-8) { std::cout << "input:" << std::endl << input << std::endl; std::cout << "output:" << std::endl << output << std::endl; std::cout << "my_output:" << std::endl << my_output << std::endl; } } BOOST_AUTO_TEST_CASE(log_m_test) { using Mat_t = Eigen::Matrix; Mat_t input{}; constexpr Real log_tol{1e-8}; input << 1.32460909, 0.86867096, 0, 0.86867096, 1.32460909, 0, 0, 0, 2.71828183; Mat_t output{}; output << 0, .25 * pi, 0, .25 * pi, 0, 0, 0, 0, 1; auto my_output{logm(input)}; Real error{(my_output - output).norm() / output.norm()}; BOOST_CHECK_LT(error, log_tol); if (error >= log_tol) { std::cout << "input:" << std::endl << input << std::endl; std::cout << "output:" << std::endl << output << std::endl; std::cout << "my_output:" << std::endl << my_output << std::endl; } input << 1.0001000000000002, 0.010000000000000116, 0, 0.010000000000000061, 1.0000000000000002, 0, 0, 0, 1; // from scipy.linalg.logm output << 4.99991667e-05, 9.99983334e-03, 0.00000000e+00, 9.99983334e-03, -4.99991667e-05, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00; my_output = logm(input); error = (my_output - output).norm() / output.norm(); BOOST_CHECK_LT(error, log_tol); if (error >= log_tol) { std::cout << "input:" << std::endl << input << std::endl; std::cout << "output:" << std::endl << output << std::endl; std::cout << "my_output:" << std::endl << my_output << std::endl; } input << 1.0001000000000002, 0.010000000000000116, 0, 0.010000000000000061, 1.0000000000000002, 0, 0, 0, 1; input = input.transpose().eval(); // from scipy.linalg.logm output << 4.99991667e-05, 9.99983334e-03, 0.00000000e+00, 9.99983334e-03, -4.99991667e-05, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00; my_output = logm(input); error = (my_output - output).norm() / output.norm(); BOOST_CHECK_LT(error, log_tol); if (error >= log_tol) { std::cout << "input:" << std::endl << input << std::endl; std::cout << "output:" << std::endl << output << std::endl; std::cout << "my_output:" << std::endl << my_output << std::endl; } Mat_t my_output_alt{logm_alt(input)}; error = (my_output_alt - output).norm() / output.norm(); BOOST_CHECK_LT(error, log_tol); if (error >= log_tol) { std::cout << "input:" << std::endl << input << std::endl; std::cout << "output:" << std::endl << output << std::endl; std::cout << "my_output:" << std::endl << my_output_alt << std::endl; } } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muGrid diff --git a/tests/header_test_fields.cc b/tests/header_test_fields.cc index 509be20..c105424 100644 --- a/tests/header_test_fields.cc +++ b/tests/header_test_fields.cc @@ -1,273 +1,273 @@ /** * @file header_test_fields.cc * * @author Till Junge * * @date 20 Sep 2017 * * @brief Test Fields that are used in FieldCollections * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "tests.hh" -#include "common/field_collection.hh" -#include "common/field.hh" -#include "common/ccoord_operations.hh" +#include +#include +#include #include #include -namespace muSpectre { +namespace muGrid { BOOST_AUTO_TEST_SUITE(field_test); template struct FieldFixture { constexpr static bool IsGlobal{Global}; constexpr static Dim_t Order{secondOrder}; constexpr static Dim_t SDim{twoD}; constexpr static Dim_t MDim{threeD}; constexpr static Dim_t NbComponents{ipow(MDim, Order)}; using FieldColl_t = std::conditional_t, LocalFieldCollection>; using TField_t = TensorField; using MField_t = MatrixField; using DField_t = TypedField; FieldFixture() : tensor_field{make_field("TensorField", this->fc)}, matrix_field{make_field("MatrixField", this->fc)}, dynamic_field1{make_field( "Dynamically sized field with correct number of" " components", this->fc, ipow(MDim, Order))}, dynamic_field2{make_field( "Dynamically sized field with incorrect number" " of components", this->fc, NbComponents + 1)} {} ~FieldFixture() = default; FieldColl_t fc{}; TField_t & tensor_field; MField_t & matrix_field; DField_t & dynamic_field1; DField_t & dynamic_field2; }; using field_fixtures = boost::mpl::list, FieldFixture>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE(size_check_global, FieldFixture) { // check that fields are initialised with empty vector BOOST_CHECK_EQUAL(tensor_field.size(), 0); BOOST_CHECK_EQUAL(dynamic_field1.size(), 0); BOOST_CHECK_EQUAL(dynamic_field2.size(), 0); // check that returned size is correct Dim_t len{2}; auto sizes{CcoordOps::get_cube(len)}; fc.initialise(sizes, {}); auto nb_pixels{CcoordOps::get_size(sizes)}; BOOST_CHECK_EQUAL(tensor_field.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field1.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field2.size(), nb_pixels); constexpr Dim_t pad_size{3}; tensor_field.set_pad_size(pad_size); dynamic_field1.set_pad_size(pad_size); dynamic_field2.set_pad_size(pad_size); // check that setting pad size won't change logical size BOOST_CHECK_EQUAL(tensor_field.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field1.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field2.size(), nb_pixels); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE(size_check_local, FieldFixture) { // check that fields are initialised with empty vector BOOST_CHECK_EQUAL(tensor_field.size(), 0); BOOST_CHECK_EQUAL(dynamic_field1.size(), 0); BOOST_CHECK_EQUAL(dynamic_field2.size(), 0); // check that returned size is correct Dim_t nb_pixels{3}; Eigen::Array new_elem; Eigen::Array wrong_elem; for (Dim_t i{0}; i < NbComponents; ++i) { new_elem(i) = i; wrong_elem(i) = .1 * i; } for (Dim_t i{0}; i < nb_pixels; ++i) { tensor_field.push_back(new_elem); dynamic_field1.push_back(new_elem); BOOST_CHECK_THROW(dynamic_field1.push_back(wrong_elem), FieldError); BOOST_CHECK_THROW(dynamic_field2.push_back(new_elem), FieldError); } BOOST_CHECK_EQUAL(tensor_field.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field1.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field2.size(), 0); for (Dim_t i{0}; i < nb_pixels; ++i) { fc.add_pixel({i}); } fc.initialise(); BOOST_CHECK_EQUAL(tensor_field.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field1.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field2.size(), nb_pixels); BOOST_CHECK_EQUAL(tensor_field.get_pad_size(), 0); BOOST_CHECK_EQUAL(dynamic_field1.get_pad_size(), 0); BOOST_CHECK_EQUAL(dynamic_field2.get_pad_size(), 0); constexpr Dim_t pad_size{3}; tensor_field.set_pad_size(pad_size); dynamic_field1.set_pad_size(pad_size); dynamic_field2.set_pad_size(pad_size); BOOST_CHECK_EQUAL(tensor_field.get_pad_size(), pad_size); BOOST_CHECK_EQUAL(dynamic_field1.get_pad_size(), pad_size); BOOST_CHECK_EQUAL(dynamic_field2.get_pad_size(), pad_size); // check that setting pad size won't change logical size BOOST_CHECK_EQUAL(tensor_field.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field1.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field2.size(), nb_pixels); } BOOST_AUTO_TEST_CASE(simple_creation) { constexpr Dim_t sdim{twoD}; constexpr Dim_t mdim{twoD}; constexpr Dim_t order{fourthOrder}; using FC_t = GlobalFieldCollection; FC_t fc; using TF_t = TensorField; auto & field{make_field("TensorField 1", fc)}; // check that fields are initialised with empty vector BOOST_CHECK_EQUAL(field.size(), 0); Dim_t len{2}; fc.initialise(CcoordOps::get_cube(len), {}); // check that returned size is correct BOOST_CHECK_EQUAL(field.size(), ipow(len, sdim)); // check that setting pad size won't change logical size field.set_pad_size(24); BOOST_CHECK_EQUAL(field.size(), ipow(len, sdim)); } BOOST_AUTO_TEST_CASE(dynamic_field_creation) { constexpr Dim_t sdim{threeD}; Dim_t nb_components{2}; using FC_t = GlobalFieldCollection; FC_t fc{}; make_field>("Dynamic Field", fc, nb_components); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(get_zeros_like, Fix, field_fixtures, Fix) { auto & t_clone{Fix::tensor_field.get_zeros_like("tensor clone")}; static_assert(std::is_same, typename Fix::TField_t>::value, "wrong overloaded function"); auto & m_clone{Fix::matrix_field.get_zeros_like("matrix clone")}; static_assert(std::is_same, typename Fix::MField_t>::value, "wrong overloaded function"); using FieldColl_t = typename Fix::FieldColl_t; using T = typename Fix::TField_t::Scalar; TypedField & t_ref{t_clone}; auto & typed_clone{t_ref.get_zeros_like("dynamically sized clone")}; static_assert(std::is_same, TypedField>::value, "Field type incorrectly deduced"); BOOST_CHECK_EQUAL(typed_clone.get_nb_components(), t_clone.get_nb_components()); auto & dyn_clone{Fix::dynamic_field1.get_zeros_like("dynamic clone")}; static_assert( std::is_same::value, "mismatch"); BOOST_CHECK_EQUAL(typed_clone.get_nb_components(), dyn_clone.get_nb_components()); } /* ---------------------------------------------------------------------- */ BOOST_AUTO_TEST_CASE(fill_global_local) { FieldFixture global; FieldFixture local; constexpr Dim_t len{2}; constexpr auto sizes{CcoordOps::get_cube::SDim>(len)}; global.fc.initialise(sizes, {}); local.fc.add_pixel({1, 1}); local.fc.add_pixel({0, 1}); local.fc.initialise(); // fill the local matrix field and then transfer it to the global field for (auto mat : local.matrix_field.get_map()) { mat.setRandom(); } global.matrix_field.fill_from_local(local.matrix_field); for (const auto & ccoord : local.fc) { const auto & a{local.matrix_field.get_map()[ccoord]}; const auto & b{global.matrix_field.get_map()[ccoord]}; const Real error{(a - b).norm()}; BOOST_CHECK_EQUAL(error, 0.); } // fill the global tensor field and then transfer it to the global field for (auto mat : global.tensor_field.get_map()) { mat.setRandom(); } local.tensor_field.fill_from_global(global.tensor_field); for (const auto & ccoord : local.fc) { const auto & a{local.matrix_field.get_map()[ccoord]}; const auto & b{global.matrix_field.get_map()[ccoord]}; const Real error{(a - b).norm()}; BOOST_CHECK_EQUAL(error, 0.); } BOOST_CHECK_THROW(local.tensor_field.fill_from_global(global.matrix_field), std::runtime_error); BOOST_CHECK_THROW(global.tensor_field.fill_from_local(local.matrix_field), std::runtime_error); } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muGrid diff --git a/tests/header_test_raw_field_map.cc b/tests/header_test_raw_field_map.cc index 40b01fb..1b079fe 100644 --- a/tests/header_test_raw_field_map.cc +++ b/tests/header_test_raw_field_map.cc @@ -1,100 +1,100 @@ /** * file header_test_raw_field_map.cc * * @author Till Junge * * @date 17 Apr 2018 * * @brief tests for the raw field map type * * @section LICENSE * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "test_field_collections.hh" -#include "common/field_map.hh" +#include -namespace muSpectre { +namespace muGrid { BOOST_AUTO_TEST_SUITE(raw_field_map_tests); BOOST_FIXTURE_TEST_CASE_TEMPLATE(iter_field_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using MSqMap = MatrixFieldMap; MSqMap Mmap{F::fc["Tensorfield Real o2"]}; auto m_it = Mmap.begin(); auto m_it_end = Mmap.end(); RawFieldMap< Eigen::Map>> raw_map{Mmap.get_field().eigenvec()}; for (auto && mat : Mmap) { mat.setRandom(); } for (auto tup : akantu::zip(Mmap, raw_map)) { auto & mat_A = std::get<0>(tup); auto & mat_B = std::get<1>(tup); BOOST_CHECK_EQUAL((mat_A - mat_B).norm(), 0.); } Mmap.get_field().eigenvec().setZero(); for (auto && mat : raw_map) { mat.setIdentity(); } for (auto && mat : Mmap) { BOOST_CHECK_EQUAL((mat - mat.Identity()).norm(), 0.); } } BOOST_AUTO_TEST_CASE(Const_correctness_test) { Eigen::VectorXd vec1(12); vec1.setRandom(); RawFieldMap> map1{vec1}; static_assert(not map1.IsConst, "should not have been const"); RawFieldMap> cmap1{vec1}; static_assert(cmap1.IsConst, "should have been const"); const Eigen::VectorXd vec2{vec1}; RawFieldMap> cmap2{vec2}; } BOOST_AUTO_TEST_CASE(incompatible_size_check) { Eigen::VectorXd vec1(11); using RawFieldMap_t = RawFieldMap>; BOOST_CHECK_THROW(RawFieldMap_t{vec1}, std::runtime_error); } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muGrid diff --git a/tests/header_test_ref_array.cc b/tests/header_test_ref_array.cc index b11e4dd..5ac6a69 100644 --- a/tests/header_test_ref_array.cc +++ b/tests/header_test_ref_array.cc @@ -1,52 +1,53 @@ /** * @file header_test_ref_array.cc * * @author Till Junge * * @date 04 Dec 2018 * * @brief tests for the RefArray convenience struct * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ -#include "common/ref_array.hh" #include "tests.hh" -namespace muSpectre { +#include + +namespace muGrid { BOOST_AUTO_TEST_SUITE(RefArray_tests); BOOST_AUTO_TEST_CASE(two_d_test) { std::array values{2, 3}; RefArray refs{values[0], values[1]}; BOOST_CHECK_EQUAL(values[0], refs[0]); BOOST_CHECK_EQUAL(values[1], refs[1]); } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muGrid diff --git a/tests/header_test_statefields.cc b/tests/header_test_statefields.cc index e409c2b..7dd87c8 100644 --- a/tests/header_test_statefields.cc +++ b/tests/header_test_statefields.cc @@ -1,299 +1,301 @@ /** * file header_test_statefields.cc * * @author Till Junge * * @date 01 Mar 2018 * * @brief Test the StateField abstraction and the associated maps * * @section LICENSE * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ -#include "common/field.hh" -#include "common/field_collection.hh" -#include "common/statefield.hh" -#include "common/ccoord_operations.hh" +#include +#include +#include +#include #include "tests.hh" #include #include +#include -namespace muSpectre { +namespace muGrid { + using muSpectre::tol; template struct SF_Fixture { using FC_t = std::conditional_t, LocalFieldCollection>; using Field_t = TensorField; using ScalField_t = ScalarField; constexpr static size_t nb_mem{2}; constexpr static Dim_t sdim{DimS}; constexpr static Dim_t mdim{DimM}; constexpr static bool global{Global}; constexpr static size_t get_nb_mem() { return nb_mem; } constexpr static Dim_t get_sdim() { return sdim; } constexpr static Dim_t get_mdim() { return mdim; } constexpr static bool get_global() { return global; } SF_Fixture() : fc{}, sf{make_statefield>("prefix", fc)}, scalar_f{ make_statefield>("scalar", fc)}, self{*this} {} FC_t fc; StateField & sf; StateField & scalar_f; SF_Fixture & self; }; using typelist = boost::mpl::list< SF_Fixture, SF_Fixture, SF_Fixture, SF_Fixture, SF_Fixture, SF_Fixture>; BOOST_AUTO_TEST_SUITE(statefield); BOOST_AUTO_TEST_CASE(old_values_test) { constexpr Dim_t Dim{twoD}; constexpr size_t NbMem{2}; constexpr bool verbose{false}; using FC_t = LocalFieldCollection; FC_t fc{}; using Field_t = ScalarField; auto & statefield{make_statefield>("name", fc)}; fc.add_pixel({}); fc.initialise(); for (size_t i{0}; i < NbMem + 1; ++i) { statefield.current().eigen() = i + 1; if (verbose) { std::cout << "current = " << statefield.current().eigen() << std::endl << "old 1 = " << statefield.old().eigen() << std::endl << "old 2 = " << statefield.template old<2>().eigen() << std::endl << "indices = " << statefield.get_indices() << std::endl << std::endl; } statefield.cycle(); } BOOST_CHECK_EQUAL(statefield.current().eigen()(0), 1); BOOST_CHECK_EQUAL(statefield.old().eigen()(0), 3); BOOST_CHECK_EQUAL(statefield.template old<2>().eigen()(0), 2); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, Fix, typelist, Fix) { const std::string ref{"prefix"}; const std::string & fix{Fix::sf.get_prefix()}; BOOST_CHECK_EQUAL(ref, fix); } namespace internal { template struct init { static void run(Fixture_t & fix) { constexpr Dim_t dim{std::remove_reference_t::sdim}; fix.fc.initialise(CcoordOps::get_cube(3), CcoordOps::get_cube(0)); } }; template struct init { static void run(Fixture_t & fix) { constexpr Dim_t dim{std::remove_reference_t::sdim}; CcoordOps::Pixels pixels(CcoordOps::get_cube(3), CcoordOps::get_cube(0)); for (auto && pix : pixels) { fix.fc.add_pixel(pix); } fix.fc.initialise(); } }; } // namespace internal BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_iteration, Fix, typelist, Fix) { internal::init::run(Fix::self); constexpr Dim_t mdim{Fix::mdim}; constexpr bool verbose{false}; using StateFMap = StateFieldMap, Fix::nb_mem>; StateFMap matrix_map(Fix::sf); for (size_t i = 0; i < Fix::nb_mem + 1; ++i) { for (auto && wrapper : matrix_map) { wrapper.current() += (i + 1) * wrapper.current().Identity(); if (verbose) { std::cout << "pixel " << wrapper.get_ccoord() << ", memory cycle " << i << std::endl; std::cout << wrapper.current() << std::endl; std::cout << wrapper.old() << std::endl; std::cout << wrapper.template old<2>() << std::endl << std::endl; } } Fix::sf.cycle(); } for (auto && wrapper : matrix_map) { auto I{wrapper.current().Identity()}; Real error{(wrapper.current() - I).norm()}; BOOST_CHECK_LT(error, tol); error = (wrapper.old() - 3 * I).norm(); BOOST_CHECK_LT(error, tol); error = (wrapper.template old<2>() - 2 * I).norm(); BOOST_CHECK_LT(error, tol); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_default_map, Fix, typelist, Fix) { internal::init::run(Fix::self); constexpr bool verbose{false}; auto matrix_map{Fix::sf.get_map()}; for (size_t i = 0; i < Fix::nb_mem + 1; ++i) { for (auto && wrapper : matrix_map) { wrapper.current() += (i + 1) * wrapper.current().Identity(); if (verbose) { std::cout << "pixel " << wrapper.get_ccoord() << ", memory cycle " << i << std::endl; std::cout << wrapper.current() << std::endl; std::cout << wrapper.old() << std::endl; std::cout << wrapper.template old<2>() << std::endl << std::endl; } } Fix::sf.cycle(); } auto matrix_const_map{Fix::sf.get_const_map()}; for (auto && wrapper : matrix_const_map) { auto I{wrapper.current().Identity()}; Real error{(wrapper.current() - I).norm()}; BOOST_CHECK_LT(error, tol); error = (wrapper.old() - 3 * I).norm(); BOOST_CHECK_LT(error, tol); error = (wrapper.template old<2>() - 2 * I).norm(); BOOST_CHECK_LT(error, tol); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_scalar_map, Fix, typelist, Fix) { internal::init::run(Fix::self); constexpr bool verbose{false}; auto scalar_map{Fix::scalar_f.get_map()}; for (size_t i = 0; i < Fix::nb_mem + 1; ++i) { for (auto && wrapper : scalar_map) { wrapper.current() += (i + 1); if (verbose) { std::cout << "pixel " << wrapper.get_ccoord() << ", memory cycle " << i << std::endl; std::cout << wrapper.current() << std::endl; std::cout << wrapper.old() << std::endl; std::cout << wrapper.template old<2>() << std::endl << std::endl; } } Fix::scalar_f.cycle(); } auto scalar_const_map{Fix::scalar_f.get_const_map()}; BOOST_CHECK_EQUAL(scalar_const_map[0].current(), scalar_const_map[1].current()); for (auto wrapper : scalar_const_map) { Real error{wrapper.current() - 1}; BOOST_CHECK_LT(error, tol); error = wrapper.old() - 3; BOOST_CHECK_LT(error, tol); error = wrapper.template old<2>() - 2; BOOST_CHECK_LT(error, tol); } } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Polymorphic_access_by_name, Fix, typelist, Fix) { internal::init::run(Fix::self); // constexpr bool verbose{true}; auto & tensor_field = Fix::fc.get_statefield("prefix"); BOOST_CHECK_EQUAL(tensor_field.get_nb_memory(), Fix::get_nb_mem()); auto & field = Fix::fc.template get_current("prefix"); BOOST_CHECK_EQUAL(field.get_nb_components(), ipow(Fix::get_mdim(), secondOrder)); BOOST_CHECK_THROW(Fix::fc.template get_current("prefix"), std::runtime_error); auto & old_field = Fix::fc.template get_old("prefix"); BOOST_CHECK_EQUAL(old_field.get_nb_components(), field.get_nb_components()); BOOST_CHECK_THROW( Fix::fc.template get_old("prefix", Fix::get_nb_mem() + 1), std::out_of_range); auto & statefield{Fix::fc.get_statefield("prefix")}; auto & typed_statefield{ Fix::fc.template get_typed_statefield("prefix")}; auto map{ArrayFieldMap( typed_statefield.get_current_field())}; for (auto arr : map) { arr.setConstant(1); } Eigen::ArrayXXd field_copy{field.eigen()}; statefield.cycle(); auto & alt_old_field{typed_statefield.get_old_field()}; Real err{(field_copy - alt_old_field.eigen()).matrix().norm() / field_copy.matrix().norm()}; BOOST_CHECK_LT(err, tol); if (not(err < tol)) { std::cout << field_copy << std::endl << std::endl << typed_statefield.get_current_field().eigen() << std::endl << std::endl << typed_statefield.get_old_field(1).eigen() << std::endl << std::endl << typed_statefield.get_old_field(2).eigen() << std::endl; } } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muGrid diff --git a/tests/header_test_t4_map.cc b/tests/header_test_t4_map.cc index 3e9a12e..6c4f244 100644 --- a/tests/header_test_t4_map.cc +++ b/tests/header_test_t4_map.cc @@ -1,200 +1,199 @@ /** * @file header_test_t4_map.cc * * @author Till Junge * * @date 20 Nov 2017 * * @brief Test the fourth-order map on second-order tensor implementation * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include #include #include - +#include #include "tests.hh" #include "test_goodies.hh" -#include namespace muGrid { BOOST_AUTO_TEST_SUITE(T4map_tests); /** * Test fixture for construction of T4Map for the time being, symmetry is not * exploited */ template struct T4_fixture { T4_fixture() : matrix{}, tensor(matrix.data()) {} EIGEN_MAKE_ALIGNED_OPERATOR_NEW; using M4 = T4Mat; using T4 = T4MatMap; constexpr static Dim_t dim{Dim}; M4 matrix; T4 tensor; }; using fix_collection = boost::mpl::list, T4_fixture>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(Simple_construction_test, F, fix_collection, F) { BOOST_CHECK_EQUAL(F::tensor.cols(), F::dim * F::dim); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(write_access_test, F, fix_collection, F) { auto & t4 = F::tensor; constexpr Dim_t dim{F::dim}; Eigen::TensorFixedSize> t4c; Eigen::Map t4c_map(t4c.data()); for (Dim_t i = 0; i < F::dim; ++i) { for (Dim_t j = 0; j < F::dim; ++j) { for (Dim_t k = 0; k < F::dim; ++k) { for (Dim_t l = 0; l < F::dim; ++l) { get(t4, i, j, k, l) = 1000 * (i + 1) + 100 * (j + 1) + 10 * (k + 1) + l + 1; t4c(i, j, k, l) = 1000 * (i + 1) + 100 * (j + 1) + 10 * (k + 1) + l + 1; } } } } for (Dim_t i = 0; i < ipow(dim, 4); ++i) { BOOST_CHECK_EQUAL(F::matrix.data()[i], t4c.data()[i]); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(numpy_compatibility, F, fix_collection, F) { auto & t4 = F::tensor; typename F::M4 numpy_ref{}; if (F::dim == twoD) { numpy_ref << 1111., 1112., 1121., 1122., 1211., 1212., 1221., 1222., 2111., 2112., 2121., 2122., 2211., 2212., 2221., 2222.; } else { numpy_ref << 1111., 1112., 1113., 1121., 1122., 1123., 1131., 1132., 1133., 1211., 1212., 1213., 1221., 1222., 1223., 1231., 1232., 1233., 1311., 1312., 1313., 1321., 1322., 1323., 1331., 1332., 1333., 2111., 2112., 2113., 2121., 2122., 2123., 2131., 2132., 2133., 2211., 2212., 2213., 2221., 2222., 2223., 2231., 2232., 2233., 2311., 2312., 2313., 2321., 2322., 2323., 2331., 2332., 2333., 3111., 3112., 3113., 3121., 3122., 3123., 3131., 3132., 3133., 3211., 3212., 3213., 3221., 3222., 3223., 3231., 3232., 3233., 3311., 3312., 3313., 3321., 3322., 3323., 3331., 3332., 3333.; } for (Dim_t i = 0; i < F::dim; ++i) { for (Dim_t j = 0; j < F::dim; ++j) { for (Dim_t k = 0; k < F::dim; ++k) { for (Dim_t l = 0; l < F::dim; ++l) { get(t4, i, j, k, l) = 1000 * (i + 1) + 100 * (j + 1) + 10 * (k + 1) + l + 1; } } } } Real error{(t4 - testGoodies::from_numpy(numpy_ref)).norm()}; BOOST_CHECK_EQUAL(error, 0); if (error != 0) { std::cout << "T4:" << std::endl << t4 << std::endl; std::cout << "reshuffled np:" << std::endl << testGoodies::from_numpy(numpy_ref) << std::endl; std::cout << "original np:" << std::endl << numpy_ref << std::endl; } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(numpy_right_trans, F, fix_collection, F) { auto & t4 = F::tensor; typename F::M4 numpy_ref{}; if (F::dim == twoD) { numpy_ref << 1111., 1121., 1112., 1122., 1211., 1221., 1212., 1222., 2111., 2121., 2112., 2122., 2211., 2221., 2212., 2222.; } else { numpy_ref << 1111., 1121., 1131., 1112., 1122., 1132., 1113., 1123., 1133., 1211., 1221., 1231., 1212., 1222., 1232., 1213., 1223., 1233., 1311., 1321., 1331., 1312., 1322., 1332., 1313., 1323., 1333., 2111., 2121., 2131., 2112., 2122., 2132., 2113., 2123., 2133., 2211., 2221., 2231., 2212., 2222., 2232., 2213., 2223., 2233., 2311., 2321., 2331., 2312., 2322., 2332., 2313., 2323., 2333., 3111., 3121., 3131., 3112., 3122., 3132., 3113., 3123., 3133., 3211., 3221., 3231., 3212., 3222., 3232., 3213., 3223., 3233., 3311., 3321., 3331., 3312., 3322., 3332., 3313., 3323., 3333.; } for (Dim_t i = 0; i < F::dim; ++i) { for (Dim_t j = 0; j < F::dim; ++j) { for (Dim_t k = 0; k < F::dim; ++k) { for (Dim_t l = 0; l < F::dim; ++l) { get(t4, i, j, k, l) = 1000 * (i + 1) + 100 * (j + 1) + 10 * (k + 1) + l + 1; } } } } Real error{(t4 - testGoodies::right_transpose(numpy_ref)).norm()}; BOOST_CHECK_EQUAL(error, 0); if (error != 0) { std::cout << "T4:" << std::endl << t4 << std::endl; std::cout << "reshuffled np:" << std::endl << testGoodies::from_numpy(numpy_ref) << std::endl; std::cout << "original np:" << std::endl << numpy_ref << std::endl; } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(assign_matrix_test, F, fix_collection, F) { decltype(F::matrix) matrix; matrix.setRandom(); F::tensor = matrix; for (Dim_t i = 0; i < ipow(F::dim, 4); ++i) { BOOST_CHECK_EQUAL(F::matrix.data()[i], matrix.data()[i]); } } BOOST_AUTO_TEST_CASE(Return_ref_from_const_test) { constexpr Dim_t dim{2}; using T = int; using M4 = Eigen::Matrix; using M4c = const Eigen::Matrix; using T4 = T4MatMap; using T4c = T4MatMap; M4 mat; mat.setRandom(); M4c cmat{mat}; T4 tensor{mat.data()}; T4c ctensor{mat.data()}; T a = get(tensor, 0, 0, 0, 1); T b = get(ctensor, 0, 0, 0, 1); BOOST_CHECK_EQUAL(a, b); } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muGrid diff --git a/tests/header_test_tensor_algebra.cc b/tests/header_test_tensor_algebra.cc index 7db7f91..6bd9998 100644 --- a/tests/header_test_tensor_algebra.cc +++ b/tests/header_test_tensor_algebra.cc @@ -1,296 +1,297 @@ /** * @file header_test_tensor_algebra.cc * * @author Till Junge * * @date 05 Nov 2017 * * @brief Tests for the tensor algebra functions * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include -#include "common/tensor_algebra.hh" #include "tests.hh" #include "tests/test_goodies.hh" +#include +#include -namespace muSpectre { - +namespace muGrid { + using muSpectre::tol; BOOST_AUTO_TEST_SUITE(tensor_algebra) auto TerrNorm = [](auto && t) { return Eigen::Tensor(t.abs().sum())(); }; /* ---------------------------------------------------------------------- */ BOOST_AUTO_TEST_CASE(tensor_outer_product_test) { constexpr Dim_t dim{2}; Eigen::TensorFixedSize> A, B; // use prime numbers so that every multiple is uniquely identifiable A.setValues({{1, 2}, {3, 7}}); B.setValues({{11, 13}, {17, 19}}); Eigen::TensorFixedSize> Res1, Res2, Res3; for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t k = 0; k < dim; ++k) { for (Dim_t l = 0; l < dim; ++l) { Res1(i, j, k, l) = A(i, j) * B(k, l); Res2(i, j, k, l) = A(i, k) * B(j, l); Res3(i, j, k, l) = A(i, l) * B(j, k); } } } } Real error = TerrNorm(Res1 - Tensors::outer(A, B)); BOOST_CHECK_LT(error, tol); error = TerrNorm(Res2 - Tensors::outer_under(A, B)); BOOST_CHECK_LT(error, tol); error = TerrNorm(Res3 - Tensors::outer_over(A, B)); if (error > tol) { std::cout << "reference:" << std::endl << Res3 << std::endl; std::cout << "result:" << std::endl << Tensors::outer_over(A, B) << std::endl; std::cout << "A:" << std::endl << A << std::endl; std::cout << "B" << std::endl << B << std::endl; decltype(Res3) tmp = Tensors::outer_over(A, B); for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t k = 0; k < dim; ++k) { for (Dim_t l = 0; l < dim; ++l) { std::cout << "for (" << i << ", " << j << ", " << k << ", " << l << "), ref: " << std::setw(3) << Res3(i, j, k, l) << ", res: " << std::setw(3) << tmp(i, j, k, l) << std::endl; } } } } } BOOST_CHECK_LT(error, tol); error = TerrNorm(Res3 - Tensors::outer(A, B)); BOOST_CHECK_GT(error, tol); }; BOOST_FIXTURE_TEST_CASE_TEMPLATE(outer_products, Fix, testGoodies::dimlist, Fix) { constexpr auto dim{Fix::dim}; using T2 = Tensors::Tens2_t; using M2 = Matrices::Tens2_t; using Map2 = Eigen::Map; using T4 = Tensors::Tens4_t; using M4 = Matrices::Tens4_t; using Map4 = Eigen::Map; T2 A, B; T4 RT; A.setRandom(); B.setRandom(); Map2 Amap(A.data()); Map2 Bmap(B.data()); M2 C, D; M4 RM; C = Amap; D = Bmap; auto error = [](const T4 & A, const M4 & B) { return (B - Map4(A.data())).norm(); }; // Check outer product RT = Tensors::outer(A, B); RM = Matrices::outer(C, D); BOOST_CHECK_LT(error(RT, RM), tol); // Check outer_under product RT = Tensors::outer_under(A, B); RM = Matrices::outer_under(C, D); BOOST_CHECK_LT(error(RT, RM), tol); // Check outer_over product RT = Tensors::outer_over(A, B); RM = Matrices::outer_over(C, D); BOOST_CHECK_LT(error(RT, RM), tol); } BOOST_AUTO_TEST_CASE(tensor_multiplication) { constexpr Dim_t dim{2}; using Strain_t = Eigen::TensorFixedSize>; Strain_t A, B; A.setValues({{1, 2}, {3, 7}}); B.setValues({{11, 13}, {17, 19}}); Strain_t FF1 = A * B; // element-wise multiplication std::array, 1> prod_dims{Eigen::IndexPair{1, 0}}; Strain_t FF2 = A.contract(B, prod_dims); // correct option 1 Strain_t FF3; using Mat_t = Eigen::Map>; // following only works for evaluated tensors (which already have data()) Mat_t(FF3.data()) = Mat_t(A.data()) * Mat_t(B.data()); Strain_t ref; ref.setZero(); for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t a = 0; a < dim; ++a) { ref(i, j) += A(i, a) * B(a, j); } } } using Strain_tw = Eigen::TensorFixedSize>; Strain_tw C; C.setConstant(100); // static_assert(!std::is_convertible::value, // "Tensors not size-protected"); if (std::is_convertible::value) { // std::cout << "this is not good, should I abandon Tensors?"; } // this test seems useless. I use to detect if Eigen changed the // default tensor product Real error = TerrNorm(FF1 - ref); if (error < tol) { std::cout << "A =" << std::endl << A << std::endl; std::cout << "B =" << std::endl << B << std::endl; std::cout << "FF1 =" << std::endl << FF1 << std::endl; std::cout << "ref =" << std::endl << ref << std::endl; } BOOST_CHECK_GT(error, tol); error = TerrNorm(FF2 - ref); if (error > tol) { std::cout << "A =" << std::endl << A << std::endl; std::cout << "B =" << std::endl << B << std::endl; std::cout << "FF2 =" << std::endl << FF2 << std::endl; std::cout << "ref =" << std::endl << ref << std::endl; } BOOST_CHECK_LT(error, tol); error = TerrNorm(FF3 - ref); if (error > tol) { std::cout << "A =" << std::endl << A << std::endl; std::cout << "B =" << std::endl << B << std::endl; std::cout << "FF3 =" << std::endl << FF3 << std::endl; std::cout << "ref =" << std::endl << ref << std::endl; } BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensmult, Fix, testGoodies::dimlist, Fix) { using Matrices::Tens2_t; using Matrices::Tens4_t; using Matrices::tensmult; constexpr Dim_t dim{Fix::dim}; using T4 = Tens4_t; using T2 = Tens2_t; using V2 = Eigen::Matrix; T4 C; C.setRandom(); T2 E; E.setRandom(); Eigen::Map Ev(E.data()); T2 R = tensmult(C, E); auto error = (Eigen::Map(R.data()) - C * Ev).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tracer, Fix, testGoodies::dimlist, Fix) { using Matrices::Tens2_t; using Matrices::Tens4_t; using Matrices::tensmult; constexpr Dim_t dim{Fix::dim}; using T2 = Tens2_t; auto tracer = Matrices::Itrac(); T2 F; F.setRandom(); auto Ftrac = tensmult(tracer, F); auto error = (Ftrac - F.trace() * F.Identity()).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_identity, Fix, testGoodies::dimlist, Fix) { using Matrices::Tens2_t; using Matrices::Tens4_t; using Matrices::tensmult; constexpr Dim_t dim{Fix::dim}; using T2 = Tens2_t; auto ident = Matrices::Iiden(); T2 F; F.setRandom(); auto Fiden = tensmult(ident, F); auto error = (Fiden - F).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_transposer, Fix, testGoodies::dimlist, Fix) { using Matrices::Tens2_t; using Matrices::Tens4_t; using Matrices::tensmult; constexpr Dim_t dim{Fix::dim}; using T2 = Tens2_t; auto trnst = Matrices::Itrns(); T2 F; F.setRandom(); auto Ftrns = tensmult(trnst, F); auto error = (Ftrns - F.transpose()).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_symmetriser, Fix, testGoodies::dimlist, Fix) { using Matrices::Tens2_t; using Matrices::Tens4_t; using Matrices::tensmult; constexpr Dim_t dim{Fix::dim}; using T2 = Tens2_t; auto symmt = Matrices::Isymm(); T2 F; F.setRandom(); auto Fsymm = tensmult(symmt, F); auto error = (Fsymm - .5 * (F + F.transpose())).norm(); BOOST_CHECK_LT(error, tol); } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muGrid diff --git a/tests/test_cell_base.cc b/tests/test_cell_base.cc index abde920..5b0c490 100644 --- a/tests/test_cell_base.cc +++ b/tests/test_cell_base.cc @@ -1,325 +1,327 @@ /** * @file test_cell_base.cc * * @author Till Junge * * @date 14 Dec 2017 * * @brief Tests for the basic cell class * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include #include "tests.hh" -#include "common/common.hh" -#include "common/iterators.hh" -#include "common/field_map.hh" +#include +#include #include "tests/test_goodies.hh" #include "cell/cell_factory.hh" #include "materials/material_linear_elastic1.hh" namespace muSpectre { + using muGrid::make_field; + using muGrid::TypedField; BOOST_AUTO_TEST_SUITE(cell_base); template struct Sizes {}; template <> struct Sizes { constexpr static Dim_t sdim{twoD}; constexpr static Ccoord_t get_resolution() { return Ccoord_t{3, 5}; } constexpr static Rcoord_t get_lengths() { return Rcoord_t{3.4, 5.8}; } }; template <> struct Sizes { constexpr static Dim_t sdim{threeD}; constexpr static Ccoord_t get_resolution() { return Ccoord_t{3, 5, 7}; } constexpr static Rcoord_t get_lengths() { return Rcoord_t{3.4, 5.8, 6.7}; } }; template struct CellBaseFixture : CellBase { constexpr static Dim_t sdim{DimS}; constexpr static Dim_t mdim{DimM}; constexpr static Formulation formulation{form}; CellBaseFixture() : CellBase{std::move( cell_input(Sizes::get_resolution(), Sizes::get_lengths(), form))} {} }; using fixlist = boost::mpl::list< CellBaseFixture, CellBaseFixture, CellBaseFixture, CellBaseFixture>; BOOST_AUTO_TEST_CASE(manual_construction) { constexpr Dim_t dim{twoD}; Ccoord_t resolutions{3, 3}; Rcoord_t lengths{2.3, 2.7}; Formulation form{Formulation::finite_strain}; - auto fft_ptr{std::make_unique>(resolutions, dim * dim)}; + auto fft_ptr{ + std::make_unique>(resolutions, dim * dim)}; auto proj_ptr{std::make_unique>( std::move(fft_ptr), lengths)}; CellBase sys{std::move(proj_ptr)}; auto sys2{make_cell(resolutions, lengths, form)}; auto sys2b{std::move(sys2)}; BOOST_CHECK_EQUAL(sys2b.size(), sys.size()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { BOOST_CHECK_THROW(fix::check_material_coverage(), std::runtime_error); BOOST_CHECK_THROW(fix::initialise(), std::runtime_error); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(add_material_test, fix, fixlist, fix) { constexpr Dim_t dim{fix::sdim}; using Material_t = MaterialLinearElastic1; auto Material_hard = std::make_unique("hard", 210e9, .33); BOOST_CHECK_NO_THROW(fix::add_material(std::move(Material_hard))); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(simple_evaluation_test, fix, fixlist, fix) { constexpr Dim_t dim{fix::sdim}; constexpr Formulation form{fix::formulation}; using Mat_t = MaterialLinearElastic1; const Real Young{210e9}, Poisson{.33}; const Real lambda{Young * Poisson / ((1 + Poisson) * (1 - 2 * Poisson))}; const Real mu{Young / (2 * (1 + Poisson))}; auto Material_hard = std::make_unique("hard", Young, Poisson); for (auto && pixel : *this) { Material_hard->add_pixel(pixel); } fix::add_material(std::move(Material_hard)); auto & F = fix::get_strain(); auto F_map = F.get_map(); // finite strain formulation expects the deformation gradient F, // while small strain expects infinitesimal strain ε for (auto grad : F_map) { switch (form) { case Formulation::finite_strain: { grad = grad.Identity(); break; } case Formulation::small_strain: { grad = grad.Zero(); break; } default: BOOST_CHECK(false); break; } } auto res_tup{fix::evaluate_stress_tangent(F)}; auto stress{std::get<0>(res_tup).get_map()}; auto tangent{std::get<1>(res_tup).get_map()}; - auto tup = - testGoodies::objective_hooke_explicit(lambda, mu, Matrices::I2()); + auto tup = muGrid::testGoodies::objective_hooke_explicit( + lambda, mu, Matrices::I2()); auto P_ref = std::get<0>(tup); for (auto mat : stress) { Real norm = (mat - P_ref).norm(); BOOST_CHECK_EQUAL(norm, 0.); } auto tan_ref = std::get<1>(tup); for (const auto tan : tangent) { Real norm = (tan - tan_ref).norm(); BOOST_CHECK_EQUAL(norm, 0.); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(evaluation_test, fix, fixlist, fix) { constexpr Dim_t dim{fix::sdim}; using Mat_t = MaterialLinearElastic1; auto Material_hard = std::make_unique("hard", 210e9, .33); auto Material_soft = std::make_unique("soft", 70e9, .3); for (auto && cnt_pixel : akantu::enumerate(*this)) { auto counter = std::get<0>(cnt_pixel); auto && pixel = std::get<1>(cnt_pixel); if (counter < 5) { Material_hard->add_pixel(pixel); } else { Material_soft->add_pixel(pixel); } } fix::add_material(std::move(Material_hard)); fix::add_material(std::move(Material_soft)); auto & F = fix::get_strain(); fix::evaluate_stress_tangent(F); fix::evaluate_stress_tangent(F); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(evaluation_test_new_interface, fix, fixlist, fix) { constexpr Dim_t dim{fix::sdim}; using Mat_t = MaterialLinearElastic1; auto Material_hard = std::make_unique("hard", 210e9, .33); auto Material_soft = std::make_unique("soft", 70e9, .3); for (auto && cnt_pixel : akantu::enumerate(*this)) { auto counter = std::get<0>(cnt_pixel); auto && pixel = std::get<1>(cnt_pixel); if (counter < 5) { Material_hard->add_pixel(pixel); } else { Material_soft->add_pixel(pixel); } } fix::add_material(std::move(Material_hard)); fix::add_material(std::move(Material_soft)); auto F_vec = fix::get_strain_vector(); F_vec.setZero(); fix::evaluate_stress_tangent(); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_managed_fields, Fix, fixlist, Fix) { Cell & dyn_handle{*this}; CellBase & base_handle{*this}; const std::string name1{"aaa"}; constexpr size_t nb_comp{5}; auto new_dyn_array{dyn_handle.get_managed_real_array(name1, nb_comp)}; BOOST_CHECK_EQUAL(new_dyn_array.rows(), nb_comp); BOOST_CHECK_EQUAL(new_dyn_array.cols(), dyn_handle.size()); BOOST_CHECK_THROW(dyn_handle.get_managed_real_array(name1, nb_comp + 1), std::runtime_error); auto & new_field{base_handle.get_managed_real_field(name1, nb_comp)}; BOOST_CHECK_EQUAL(new_field.get_nb_components(), nb_comp); BOOST_CHECK_EQUAL(new_field.size(), dyn_handle.size()); BOOST_CHECK_THROW(base_handle.get_managed_real_field(name1, nb_comp + 1), std::runtime_error); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_globalised_fields, Fix, fixlist, Fix) { constexpr Dim_t Dim{Fix::sdim}; using Mat_t = MaterialLinearElastic1; using LColl_t = typename Mat_t::MFieldCollection_t; auto & material_soft{Mat_t::make(*this, "soft", 70e9, .3)}; auto & material_hard{Mat_t::make(*this, "hard", 210e9, .3)}; for (auto && tup : akantu::enumerate(*this)) { const auto & i{std::get<0>(tup)}; const auto & pixel{std::get<1>(tup)}; if (i % 2) { material_soft.add_pixel(pixel); } else { material_hard.add_pixel(pixel); } } material_soft.initialise(); material_hard.initialise(); auto & col_soft{material_soft.get_collection()}; auto & col_hard{material_hard.get_collection()}; // compatible fields: const std::string compatible_name{"compatible"}; auto & compatible_soft{ make_field>(compatible_name, col_soft, Dim)}; auto & compatible_hard{ make_field>(compatible_name, col_hard, Dim)}; auto pixler = [](auto & field) { auto map{field.get_map()}; for (auto && tup : map.enumerate()) { const auto & pixel{std::get<0>(tup)}; auto & val{std::get<1>(tup)}; for (Dim_t i{0}; i < Dim; ++i) { val(i) = pixel[i]; } } }; pixler(compatible_soft); pixler(compatible_hard); auto & global_compatible_field{ this->get_globalised_internal_real_field(compatible_name)}; auto glo_map{global_compatible_field.get_map()}; for (auto && tup : glo_map.enumerate()) { const auto & pixel{std::get<0>(tup)}; const auto & val(std::get<1>(tup)); using Map_t = Eigen::Map>; Real err{ (val - Map_t(pixel.data()).template cast()).matrix().norm()}; BOOST_CHECK_LT(err, tol); } // incompatible fields: const std::string incompatible_name{"incompatible"}; make_field>(incompatible_name, col_soft, Dim); make_field>(incompatible_name, col_hard, Dim + 1); BOOST_CHECK_THROW( this->get_globalised_internal_real_field(incompatible_name), std::runtime_error); // wrong name/ inexistant field const std::string wrong_name{"wrong_name"}; BOOST_CHECK_THROW(this->get_globalised_internal_real_field(wrong_name), std::runtime_error); // wrong scalar type: const std::string wrong_scalar_name{"wrong_scalar"}; make_field>(wrong_scalar_name, col_soft, Dim); make_field>(wrong_scalar_name, col_hard, Dim); BOOST_CHECK_THROW( this->get_globalised_internal_real_field(wrong_scalar_name), std::runtime_error); } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_fft_utils.cc b/tests/test_fft_utils.cc index f7e12da..a9c280a 100644 --- a/tests/test_fft_utils.cc +++ b/tests/test_fft_utils.cc @@ -1,88 +1,90 @@ /** * @file test_fft_utils.cc * * @author Till Junge * * @date 11 Dec 2017 * * @brief test the small utility functions used by the fft engines and * projections * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include "tests.hh" -#include "projection/fft_utils.hh" +#include -namespace muSpectre { +namespace muFFT { BOOST_AUTO_TEST_SUITE(fft_utils); BOOST_AUTO_TEST_CASE(fft_freqs_test) { // simply comparing to np.fft.fftfreq(12, 1/12.) const std::valarray ref{0., 1., 2., 3., 4., 5., -6., -5., -4., -3., -2., -1.}; auto res{fft_freqs(12)}; Real error = std::abs(res - ref).sum(); BOOST_CHECK_EQUAL(error, 0.); } BOOST_AUTO_TEST_CASE(fft_freqs_test_length) { // simply comparing to np.fft.fftfreq(10) const std::valarray ref{0., 0.1, 0.2, 0.3, 0.4, -0.5, -0.4, -0.3, -0.2, -0.1}; auto res{fft_freqs(10, 10.)}; Real error = std::abs(res - ref).sum(); BOOST_CHECK_EQUAL(error, 0.); } BOOST_AUTO_TEST_CASE(wave_vector_computation) { // here, build a FFT_freqs and check it returns the correct xi's constexpr Dim_t dim{twoD}; FFT_freqs freq_struc{{12, 10}, {1., 10.}}; Ccoord_t ccoord1{2, 3}; auto xi{freq_struc.get_xi(ccoord1)}; auto unit_xi{freq_struc.get_unit_xi(ccoord1)}; typename FFT_freqs::Vector ref; ref << 2., .3; // from above tests - BOOST_CHECK_LT((xi - ref).norm(), tol); - BOOST_CHECK_LT(std::abs(xi.dot(unit_xi) - xi.norm()), xi.norm() * tol); - BOOST_CHECK_LT(std::abs(unit_xi.norm() - 1.), tol); + BOOST_CHECK_LT((xi - ref).norm(), muSpectre::tol); + BOOST_CHECK_LT(std::abs(xi.dot(unit_xi) - xi.norm()), + xi.norm() * muSpectre::tol); + BOOST_CHECK_LT(std::abs(unit_xi.norm() - 1.), muSpectre::tol); ccoord1 = {7, 8}; xi = freq_struc.get_xi(ccoord1); unit_xi = freq_struc.get_unit_xi(ccoord1); ref << -5., -.2; - BOOST_CHECK_LT((xi - ref).norm(), tol); - BOOST_CHECK_LT(std::abs(xi.dot(unit_xi) - xi.norm()), xi.norm() * tol); - BOOST_CHECK_LT(std::abs(unit_xi.norm() - 1.), tol); + BOOST_CHECK_LT((xi - ref).norm(), muSpectre::tol); + BOOST_CHECK_LT(std::abs(xi.dot(unit_xi) - xi.norm()), + xi.norm() * muSpectre::tol); + BOOST_CHECK_LT(std::abs(unit_xi.norm() - 1.), muSpectre::tol); } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muFFT diff --git a/tests/test_fftw_engine.cc b/tests/test_fftw_engine.cc index 05425e2..0d5acbf 100644 --- a/tests/test_fftw_engine.cc +++ b/tests/test_fftw_engine.cc @@ -1,144 +1,151 @@ /** * @file test_fftw_engine.cc * * @author Till Junge * * @date 05 Dec 2017 * * @brief tests for the fftw fft engine implementation * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include "tests.hh" -#include "projection/fftw_engine.hh" -#include "common/ccoord_operations.hh" -#include "common/field_collection.hh" -#include "common/field_map.hh" -#include "common/iterators.hh" +#include +#include +#include +#include +#include +#include -namespace muSpectre { +namespace muFFT { + using muSpectre::tol; BOOST_AUTO_TEST_SUITE(fftw_engine); /* ---------------------------------------------------------------------- */ template struct FFTW_fixture { constexpr static Dim_t box_resolution{resolution}; constexpr static Real box_length{4.5}; constexpr static Dim_t sdim{DimS}; constexpr static Dim_t mdim{DimM}; constexpr static Ccoord_t res() { - return CcoordOps::get_cube(box_resolution); + return muGrid::CcoordOps::get_cube(box_resolution); } constexpr static Ccoord_t loc() { - return CcoordOps::get_cube(0); + return muGrid::CcoordOps::get_cube(0); } FFTW_fixture() : engine(res(), DimM * DimM) {} FFTWEngine engine; }; struct FFTW_fixture_python_segfault { constexpr static Dim_t dim{twoD}; constexpr static Dim_t sdim{twoD}; constexpr static Dim_t mdim{twoD}; constexpr static Ccoord_t res() { return {6, 4}; } constexpr static Ccoord_t loc() { return {0, 0}; } FFTW_fixture_python_segfault() : engine{res(), mdim * mdim} {} FFTWEngine engine; }; using fixlist = boost::mpl::list< FFTW_fixture, FFTW_fixture, FFTW_fixture, FFTW_fixture, FFTW_fixture, FFTW_fixture, FFTW_fixture_python_segfault>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Constructor_test, Fix, fixlist, Fix) { BOOST_CHECK_NO_THROW(Fix::engine.initialise(FFT_PlanFlags::estimate)); - BOOST_CHECK_EQUAL(Fix::engine.size(), CcoordOps::get_size(Fix::res())); + BOOST_CHECK_EQUAL(Fix::engine.size(), + muGrid::CcoordOps::get_size(Fix::res())); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(fft_test, Fix, fixlist, Fix) { Fix::engine.initialise(FFT_PlanFlags::estimate); constexpr Dim_t order{2}; - using FC_t = GlobalFieldCollection; + using FC_t = muGrid::GlobalFieldCollection; FC_t fc; auto & input{ - make_field>("input", fc)}; + muGrid::make_field>( + "input", fc)}; auto & ref{ - make_field>("reference", fc)}; + muGrid::make_field>( + "reference", fc)}; auto & result{ - make_field>("result", fc)}; + muGrid::make_field>( + "result", fc)}; fc.initialise(Fix::res(), Fix::loc()); - using map_t = MatrixFieldMap; + using map_t = muGrid::MatrixFieldMap; map_t inmap{input}; auto refmap{map_t{ref}}; auto resultmap{map_t{result}}; size_t cntr{0}; for (auto tup : akantu::zip(inmap, refmap)) { cntr++; auto & in_{std::get<0>(tup)}; auto & ref_{std::get<1>(tup)}; in_.setRandom(); ref_ = in_; } auto & complex_field = Fix::engine.fft(input); - using cmap_t = MatrixFieldMap, Complex, - Fix::mdim, Fix::mdim>; + using cmap_t = + muGrid::MatrixFieldMap, Complex, + Fix::mdim, Fix::mdim>; cmap_t complex_map(complex_field); Real error = complex_map[0].imag().norm(); BOOST_CHECK_LT(error, tol); /* make sure, the engine has not modified input (which is unfortunately const-casted internally, hence this test) */ for (auto && tup : akantu::zip(inmap, refmap)) { Real error{(std::get<0>(tup) - std::get<1>(tup)).norm()}; BOOST_CHECK_LT(error, tol); } /* make sure that the ifft of fft returns the original*/ Fix::engine.ifft(result); for (auto && tup : akantu::zip(resultmap, refmap)) { Real error{ (std::get<0>(tup) * Fix::engine.normalisation() - std::get<1>(tup)) .norm()}; BOOST_CHECK_LT(error, tol); if (error > tol) { std::cout << std::get<0>(tup).array() / std::get<1>(tup).array() << std::endl << std::endl; } } } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muFFT diff --git a/tests/test_field_collections_1.cc b/tests/test_field_collections_1.cc index cb3260c..9d4a7a3 100644 --- a/tests/test_field_collections_1.cc +++ b/tests/test_field_collections_1.cc @@ -1,222 +1,222 @@ /** * @file test_field_collections_1.cc * * @author Till Junge * * @date 20 Sep 2017 * * @brief Test the FieldCollection classes which provide fast optimized * iterators over run-time typed fields * * 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 µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "test_field_collections.hh" -namespace muSpectre { +namespace muGrid { BOOST_AUTO_TEST_SUITE(field_collection_tests); BOOST_AUTO_TEST_CASE(simple) { constexpr Dim_t sdim = 2; using FC_t = GlobalFieldCollection; FC_t fc; BOOST_CHECK_EQUAL(FC_t::spatial_dim(), sdim); BOOST_CHECK_EQUAL(fc.get_spatial_dim(), sdim); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(Simple_construction_test, F, test_collections, F) { BOOST_CHECK_EQUAL(F::FC_t::spatial_dim(), F::sdim()); BOOST_CHECK_EQUAL(F::fc.get_spatial_dim(), F::sdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(get_field2_test, F, test_collections, F) { const auto order{2}; using FC_t = typename F::FC_t; using TF_t = TensorField; auto && myfield = make_field("TensorField real 2", F::fc); using TensorMap = TensorFieldMap; using MatrixMap = MatrixFieldMap; using ArrayMap = ArrayFieldMap; TensorMap TFM(myfield); MatrixMap MFM(myfield); ArrayMap AFM(myfield); BOOST_CHECK_EQUAL(TFM.info_string(), "Tensor(d, " + std::to_string(order) + "_o, " + std::to_string(F::mdim()) + "_d)"); BOOST_CHECK_EQUAL(MFM.info_string(), "Matrix(d, " + std::to_string(F::mdim()) + "x" + std::to_string(F::mdim()) + ")"); BOOST_CHECK_EQUAL(AFM.info_string(), "Array(d, " + std::to_string(F::mdim()) + "x" + std::to_string(F::mdim()) + ")"); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(multi_field_test, F, mult_collections, F) { using FC_t = typename F::FC_t; // possible maptypes for Real tensor fields using T_type = Real; using T_TFM1_t = TensorFieldMap; using T_TFM2_t = TensorFieldMap; //! dangerous using T4_Map_t = T4MatrixFieldMap; // impossible maptypes for Real tensor fields using T_SFM_t = ScalarFieldMap; using T_MFM_t = MatrixFieldMap; using T_AFM_t = ArrayFieldMap; using T_MFMw1_t = MatrixFieldMap; using T_MFMw2_t = MatrixFieldMap; using T_MFMw3_t = MatrixFieldMap; const std::string T_name{"Tensorfield Real o4"}; const std::string T_name_w{"TensorField Real o4 wrongname"}; BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_NO_THROW(T_TFM1_t(F::fc.at(T_name))); BOOST_CHECK_NO_THROW(T_TFM2_t(F::fc.at(T_name))); BOOST_CHECK_NO_THROW(T4_Map_t(F::fc.at(T_name))); BOOST_CHECK_THROW(T4_Map_t(F::fc.at(T_name_w)), std::out_of_range); BOOST_CHECK_THROW(T_MFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_AFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw1_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw3_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name_w)), std::out_of_range); // possible maptypes for integer scalar fields using S_type = Int; using S_SFM_t = ScalarFieldMap; using S_TFM1_t = TensorFieldMap; using S_TFM2_t = TensorFieldMap; using S_MFM_t = MatrixFieldMap; using S_AFM_t = ArrayFieldMap; using S4_Map_t = T4MatrixFieldMap; // impossible maptypes for integer scalar fields using S_MFMw1_t = MatrixFieldMap; using S_MFMw2_t = MatrixFieldMap; using S_MFMw3_t = MatrixFieldMap; const std::string S_name{"integer Scalar"}; const std::string S_name_w{"integer Scalar wrongname"}; BOOST_CHECK_NO_THROW(S_SFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_TFM1_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_TFM2_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_MFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_AFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S4_Map_t(F::fc.at(S_name))); BOOST_CHECK_THROW(S_MFMw1_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(T4_Map_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw3_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_SFM_t(F::fc.at(S_name_w)), std::out_of_range); // possible maptypes for complex matrix fields using M_type = Complex; using M_MFM_t = MatrixFieldMap; using M_AFM_t = ArrayFieldMap; // impossible maptypes for complex matrix fields using M_SFM_t = ScalarFieldMap; using M_MFMw1_t = MatrixFieldMap; using M_MFMw2_t = MatrixFieldMap; using M_MFMw3_t = MatrixFieldMap; const std::string M_name{"Matrixfield Complex sdim x mdim"}; const std::string M_name_w{"Matrixfield Complex sdim x mdim wrongname"}; BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_NO_THROW(M_MFM_t(F::fc.at(M_name))); BOOST_CHECK_NO_THROW(M_AFM_t(F::fc.at(M_name))); BOOST_CHECK_THROW(M_MFMw1_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw3_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name_w)), std::out_of_range); } /* ---------------------------------------------------------------------- */ //! Check whether fields can be initialized using mult_collections_t = boost::mpl::list, FC_multi_fixture<2, 3, true>, FC_multi_fixture<3, 3, true>>; using mult_collections_f = boost::mpl::list, FC_multi_fixture<2, 3, false>, FC_multi_fixture<3, 3, false>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_glob, F, mult_collections_t, F) { Ccoord_t size; Ccoord_t loc{}; for (auto && s : size) { s = 3; } BOOST_CHECK_NO_THROW(F::fc.initialise(size, loc)); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca, F, mult_collections_f, F) { testGoodies::RandRange rng; for (int i = 0; i < 7; ++i) { Ccoord_t pixel; for (auto && s : pixel) { s = rng.randval(0, 7); } F::fc.add_pixel(pixel); } BOOST_CHECK_NO_THROW(F::fc.initialise()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca_with_push_back, F, mult_collections_f, F) { constexpr auto mdim{F::mdim()}; constexpr int nb_pix{7}; testGoodies::RandRange rng{}; using ftype = internal::TypedSizedFieldBase; using stype = Eigen::Array; auto & field = static_cast(F::fc["Tensorfield Real o4"]); field.push_back(stype()); for (int i = 0; i < nb_pix; ++i) { Ccoord_t pixel; for (auto && s : pixel) { s = rng.randval(0, 7); } F::fc.add_pixel(pixel); } BOOST_CHECK_THROW(F::fc.initialise(), FieldCollectionError); for (int i = 0; i < nb_pix - 1; ++i) { field.push_back(stype()); } BOOST_CHECK_NO_THROW(F::fc.initialise()); } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muGrid diff --git a/tests/test_field_collections_2.cc b/tests/test_field_collections_2.cc index a5645c8..2d947fe 100644 --- a/tests/test_field_collections_2.cc +++ b/tests/test_field_collections_2.cc @@ -1,308 +1,310 @@ /** * @file test_field_collections_2.cc * * @author Till Junge * * @date 23 Nov 2017 * * @brief Continuation of tests from test_field_collection.cc, split for faster * compilation * * 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 µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "test_field_collections.hh" -#include "common/field_map_dynamic.hh" +#include + +namespace muGrid { + using muSpectre::tol; -namespace muSpectre { BOOST_AUTO_TEST_SUITE(field_collection_tests); BOOST_FIXTURE_TEST_CASE_TEMPLATE(iter_field_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using Tensor4Map = TensorFieldMap; Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; TypedFieldMap dyn_map{F::fc["Tensorfield Real o4"]}; F::fc["Tensorfield Real o4"].set_zero(); for (auto && tens : T4map) { BOOST_CHECK_EQUAL(Real(Eigen::Tensor(tens.abs().sum().eval())()), 0); } for (auto && tens : T4map) { tens.setRandom(); } for (auto && tup : akantu::zip(T4map, dyn_map)) { auto & tens = std::get<0>(tup); auto & dyn = std::get<1>(tup); constexpr Dim_t nb_comp{ipow(F::mdim(), order)}; Eigen::Map> tens_arr(tens.data()); Real error{(dyn - tens_arr).matrix().norm()}; BOOST_CHECK_EQUAL(error, 0); } using Tensor2Map = TensorFieldMap; using MSqMap = MatrixFieldMap; using ASqMap = ArrayFieldMap; using A2Map = ArrayFieldMap; using WrongMap = ArrayFieldMap; Tensor2Map T2map{F::fc["Tensorfield Real o2"]}; MSqMap Mmap{F::fc["Tensorfield Real o2"]}; ASqMap Amap{F::fc["Tensorfield Real o2"]}; A2Map DynMap{F::fc["Dynamically sized Field"]}; auto & fc_ref{F::fc}; BOOST_CHECK_THROW(WrongMap{fc_ref["Dynamically sized Field"]}, FieldInterpretationError); auto t2_it = T2map.begin(); auto t2_it_end = T2map.end(); auto m_it = Mmap.begin(); auto a_it = Amap.begin(); for (; t2_it != t2_it_end; ++t2_it, ++m_it, ++a_it) { t2_it->setRandom(); auto && m = *m_it; bool comp = (m == a_it->matrix()); BOOST_CHECK(comp); } size_t counter{0}; for (auto val : DynMap) { ++counter; val += val.Ones() * counter; } counter = 0; for (auto val : DynMap) { ++counter; val -= val.Ones() * counter; auto error{val.matrix().norm()}; BOOST_CHECK_LT(error, tol); } using ScalarMap = ScalarFieldMap; ScalarMap s_map{F::fc["integer Scalar"]}; for (Uint i = 0; i < s_map.size(); ++i) { s_map[i] = i; } counter = 0; for (const auto & val : s_map) { BOOST_CHECK_EQUAL(counter++, val); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(ccoord_indexing_test, F, glob_iter_colls, F) { using FC_t = typename F::Parent::FC_t; using ScalarMap = ScalarFieldMap; ScalarMap s_map{F::fc["integer Scalar"]}; for (Uint i = 0; i < s_map.size(); ++i) { s_map[i] = i; } for (size_t i = 0; i < CcoordOps::get_size(F::fc.get_sizes()); ++i) { BOOST_CHECK_EQUAL( CcoordOps::get_index(F::fc.get_sizes(), F::fc.get_locations(), CcoordOps::get_ccoord(F::fc.get_sizes(), F::fc.get_locations(), i)), i); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(iterator_methods_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using Tensor4Map = TensorFieldMap; Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; using it_t = typename Tensor4Map::iterator; std::ptrdiff_t diff{ 3}; // arbitrary, as long as it is smaller than the container size // check constructors auto itstart = T4map.begin(); // standard way of obtaining iterator auto itend = T4map.end(); // ditto it_t it1{T4map}; it_t it2{T4map, false}; it_t it3{T4map, size_t(diff)}; BOOST_CHECK(itstart == itstart); BOOST_CHECK(itstart != itend); BOOST_CHECK_EQUAL(itstart, it1); BOOST_CHECK_EQUAL(itend, it2); // check ostream operator std::stringstream response; response << it3; BOOST_CHECK_EQUAL( response.str(), std::string("iterator on field 'Tensorfield Real o4', entry ") + std::to_string(diff)); // check move and assigment constructor (and operator+) it_t it_repl{T4map}; it_t itmove = std::move(T4map.begin()); it_t it4 = it1 + diff; it_t it7 = it4 - diff; BOOST_CHECK_EQUAL(itmove, it1); BOOST_CHECK_EQUAL(it4, it3); BOOST_CHECK_EQUAL(it7, it1); // check increments/decrements BOOST_CHECK_EQUAL(it1++, it_repl); // post-increment BOOST_CHECK_EQUAL(it1, it_repl + 1); BOOST_CHECK_EQUAL(--it1, it_repl); // pre -decrement BOOST_CHECK_EQUAL(++it1, it_repl + 1); // pre -increment BOOST_CHECK_EQUAL(it1--, it_repl + 1); // post-decrement BOOST_CHECK_EQUAL(it1, it_repl); // dereference and member-of-pointer check Eigen::Tensor Tens = *it1; Eigen::Tensor Tens2 = *itstart; Eigen::Tensor check = (Tens == Tens2).all(); BOOST_CHECK_EQUAL(static_cast(check()), true); BOOST_CHECK_NO_THROW(itstart->setZero()); // check access subscripting auto T3a = *it3; auto T3b = itstart[diff]; BOOST_CHECK( static_cast(Eigen::Tensor((T3a == T3b).all())())); // div. comparisons BOOST_CHECK_LT(itstart, itend); BOOST_CHECK(!(itend < itstart)); BOOST_CHECK(!(itstart < itstart)); BOOST_CHECK_LE(itstart, itend); BOOST_CHECK_LE(itstart, itstart); BOOST_CHECK(!(itend <= itstart)); BOOST_CHECK_GT(itend, itstart); BOOST_CHECK(!(itend > itend)); BOOST_CHECK(!(itstart > itend)); BOOST_CHECK_GE(itend, itstart); BOOST_CHECK_GE(itend, itend); BOOST_CHECK(!(itstart >= itend)); // check assignment increment/decrement BOOST_CHECK_EQUAL(it1 += diff, it3); BOOST_CHECK_EQUAL(it1 -= diff, itstart); // check cell coordinates using Ccoord = Ccoord_t; Ccoord a{itstart.get_ccoord()}; Ccoord b{0}; // Weirdly, boost::has_left_shift::value is false for // Ccoords, even though the operator is implemented :(, hence we can't write // BOOST_CHECK_EQUAL(a, b); // Workaround: bool check2 = (a == b); BOOST_CHECK(check2); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_tensor_iter_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using Tensor4Map = TensorFieldMap; Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; using T_t = typename Tensor4Map::T_t; Eigen::TensorMap Tens2(T4map[0].data(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim()); for (auto it = T4map.cbegin(); it != T4map.cend(); ++it) { // maps to const tensors can't be initialised with a const pointer this // sucks auto && tens = *it; auto && ptr = tens.data(); static_assert( std::is_pointer>::value, "should be getting a pointer"); // static_assert(std::is_const>::value, // "should be const"); // If Tensor were written well, above static_assert should pass, and the // following check shouldn't. If it get triggered, it means that a newer // version of Eigen now does have const-correct // TensorMap. This means that const-correct field maps // are then also possible for tensors BOOST_CHECK(!std::is_const>::value); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_matrix_iter_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using MatrixMap = MatrixFieldMap; MatrixMap Mmap{F::fc["Matrixfield Complex sdim x mdim"]}; for (auto it = Mmap.cbegin(); it != Mmap.cend(); ++it) { // maps to const tensors can't be initialised with a const pointer this // sucks auto && mat = *it; auto && ptr = mat.data(); static_assert( std::is_pointer>::value, "should be getting a pointer"); // static_assert(std::is_const>::value, // "should be const"); // If Matrix were written well, above static_assert should pass, and the // following check shouldn't. If it get triggered, it means that a newer // version of Eigen now does have const-correct // MatrixMap. This means that const-correct field maps // are then also possible for matrices BOOST_CHECK(!std::is_const>::value); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_scalar_iter_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using ScalarMap = ScalarFieldMap; ScalarMap Smap{F::fc["integer Scalar"]}; for (auto it = Smap.cbegin(); it != Smap.cend(); ++it) { auto && scal = *it; static_assert( std::is_const>::value, "referred type should be const"); static_assert(std::is_lvalue_reference::value, "Should have returned an lvalue ref"); } } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muGrid diff --git a/tests/test_field_collections_3.cc b/tests/test_field_collections_3.cc index 0fd7e5d..092191d 100644 --- a/tests/test_field_collections_3.cc +++ b/tests/test_field_collections_3.cc @@ -1,163 +1,163 @@ /** * @file test_field_collections_3.cc * * @author Till Junge * * @date 19 Dec 2017 * * @brief Continuation of tests from test_field_collection_2.cc, split for * faster compilation * * 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 µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "test_field_collections.hh" #include "tests/test_goodies.hh" -#include "common/tensor_algebra.hh" +#include #include -namespace muSpectre { +namespace muGrid { BOOST_AUTO_TEST_SUITE(field_collection_tests); /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(assignment_test, Fix, iter_collections, Fix) { auto t4map = Fix::t4_field.get_map(); auto t2map = Fix::t2_field.get_map(); auto scmap = Fix::sc_field.get_map(); auto m2map = Fix::m2_field.get_map(); const auto t4map_val{Matrices::Isymm()}; t4map = t4map_val; const auto t2map_val{Matrices::I2()}; t2map = t2map_val; const Int scmap_val{1}; scmap = scmap_val; Eigen::Matrix m2map_val; m2map_val.setRandom(); m2map = m2map_val; const size_t nb_pts{Fix::fc.size()}; testGoodies::RandRange rnd{}; BOOST_CHECK_EQUAL((t4map[rnd.randval(0, nb_pts - 1)] - t4map_val).norm(), 0.); BOOST_CHECK_EQUAL((t2map[rnd.randval(0, nb_pts - 1)] - t2map_val).norm(), 0.); BOOST_CHECK_EQUAL((scmap[rnd.randval(0, nb_pts - 1)] - scmap_val), 0.); BOOST_CHECK_EQUAL((m2map[rnd.randval(0, nb_pts - 1)] - m2map_val).norm(), 0.); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Eigentest, Fix, iter_collections, Fix) { auto t4eigen = Fix::t4_field.eigen(); auto t2eigen = Fix::t2_field.eigen(); BOOST_CHECK_EQUAL(t4eigen.rows(), ipow(Fix::mdim(), 4)); BOOST_CHECK_EQUAL(t4eigen.cols(), Fix::t4_field.size()); using T2_t = typename Eigen::Matrix; T2_t test_mat; test_mat.setRandom(); Eigen::Map> test_map( test_mat.data()); t2eigen.col(0) = test_map; BOOST_CHECK_EQUAL((Fix::t2_field.get_map()[0] - test_mat).norm(), 0.); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(field_proxy_test, Fix, iter_collections, Fix) { Eigen::VectorXd t4values{Fix::t4_field.eigenvec()}; using FieldProxy_t = TypedField; //! create a field proxy FieldProxy_t proxy("proxy to 'Tensorfield Real o4'", Fix::fc, t4values, Fix::t4_field.get_nb_components()); Eigen::VectorXd wrong_size_not_multiple{ Eigen::VectorXd::Zero(t4values.size() + 1)}; BOOST_CHECK_THROW(FieldProxy_t("size not a multiple of nb_components", Fix::fc, wrong_size_not_multiple, Fix::t4_field.get_nb_components()), FieldError); Eigen::VectorXd wrong_size_but_multiple{Eigen::VectorXd::Zero( t4values.size() + Fix::t4_field.get_nb_components())}; BOOST_CHECK_THROW(FieldProxy_t("size wrong multiple of nb_components", Fix::fc, wrong_size_but_multiple, Fix::t4_field.get_nb_components()), FieldError); using Tensor4Map = T4MatrixFieldMap; Tensor4Map ref_map{Fix::t4_field}; Tensor4Map proxy_map{proxy}; for (auto tup : akantu::zip(ref_map, proxy_map)) { auto & ref = std::get<0>(tup); auto & prox = std::get<1>(tup); BOOST_CHECK_EQUAL((ref - prox).norm(), 0); } } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(field_proxy_of_existing_field, Fix, iter_collections, Fix) { Eigen::Ref t4values{Fix::t4_field.eigenvec()}; using FieldProxy_t = TypedField; //! create a field proxy FieldProxy_t proxy("proxy to 'Tensorfield Real o4'", Fix::fc, t4values, Fix::t4_field.get_nb_components()); using Tensor4Map = T4MatrixFieldMap; Tensor4Map ref_map{Fix::t4_field}; Tensor4Map proxy_map{proxy}; for (auto tup : akantu::zip(ref_map, proxy_map)) { auto & ref = std::get<0>(tup); auto & prox = std::get<1>(tup); prox += prox.Identity(); BOOST_CHECK_EQUAL((ref - prox).norm(), 0); } } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(typed_field_getter, Fix, mult_collections, Fix) { constexpr auto mdim{Fix::mdim()}; auto & fc{Fix::fc}; auto & field = fc.template get_typed_field("Tensorfield Real o4"); BOOST_CHECK_EQUAL(field.get_nb_components(), ipow(mdim, fourthOrder)); } BOOST_AUTO_TEST_SUITE_END(); -} // namespace muSpectre +} // namespace muGrid diff --git a/tests/test_geometry.cc b/tests/test_geometry.cc index 1beb9cd..7f3c643 100644 --- a/tests/test_geometry.cc +++ b/tests/test_geometry.cc @@ -1,182 +1,182 @@ /** * file test_geometry.cc * * @author Till Junge * * @date 19 Apr 2018 * * @brief Tests for tensor rotations * * @section LICENSE * * 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 µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/geometry.hh" #include "tests.hh" #include "test_goodies.hh" -#include "common/T4_map_proxy.hh" +#include #include #include #include +#include namespace muSpectre { - BOOST_AUTO_TEST_SUITE(geometry); /* ---------------------------------------------------------------------- */ template struct RotationFixture { static constexpr Dim_t Dim{Dim_}; using Vec_t = Eigen::Matrix; using Mat_t = Eigen::Matrix; - using Ten_t = T4Mat; + using Ten_t = muGrid::T4Mat; using Angles_t = Eigen::Matrix; using Rot_t = Rotator; static constexpr RotationOrder EulerOrder{Rot}; static constexpr Dim_t get_Dim() { return Dim_; } RotationFixture() : rotator{euler} {} - testGoodies::RandRange rr{}; - Angles_t euler{2 * pi * Angles_t::Random()}; + muGrid::testGoodies::RandRange rr{}; + Angles_t euler{2 * muGrid::pi * Angles_t::Random()}; Vec_t v{Vec_t::Random()}; Mat_t m{Mat_t::Random()}; Ten_t t{Ten_t::Random()}; Rot_t rotator; }; /* ---------------------------------------------------------------------- */ using fix_list = boost::mpl::list, RotationFixture>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(rotation_test, Fix, fix_list, Fix) { using Vec_t = typename Fix::Vec_t; using Mat_t = typename Fix::Mat_t; using Ten_t = typename Fix::Ten_t; constexpr const Dim_t Dim{Fix::get_Dim()}; Vec_t & v{Fix::v}; Mat_t & m{Fix::m}; Ten_t & t{Fix::t}; const Mat_t & R{Fix::rotator.get_rot_mat()}; Vec_t v_ref{R * v}; Mat_t m_ref{R * m * R.transpose()}; Ten_t t_ref{Ten_t::Zero()}; for (int i = 0; i < Dim; ++i) { for (int a = 0; a < Dim; ++a) { for (int l = 0; l < Dim; ++l) { for (int b = 0; b < Dim; ++b) { for (int m = 0; m < Dim; ++m) { for (int n = 0; n < Dim; ++n) { for (int o = 0; o < Dim; ++o) { for (int p = 0; p < Dim; ++p) { - get(t_ref, a, b, o, p) += R(a, i) * R(b, l) * - get(t, i, l, m, n) * R(o, m) * - R(p, n); + muGrid::get(t_ref, a, b, o, p) += + R(a, i) * R(b, l) * muGrid::get(t, i, l, m, n) * + R(o, m) * R(p, n); } } } } } } } } Vec_t v_rotator(Fix::rotator.rotate(v)); Mat_t m_rotator(Fix::rotator.rotate(m)); Ten_t t_rotator(Fix::rotator.rotate(t)); auto v_error{(v_rotator - v_ref).norm() / v_ref.norm()}; BOOST_CHECK_LT(v_error, tol); auto m_error{(m_rotator - m_ref).norm() / m_ref.norm()}; BOOST_CHECK_LT(m_error, tol); auto t_error{(t_rotator - t_ref).norm() / t_ref.norm()}; BOOST_CHECK_LT(t_error, tol); if (t_error >= tol) { std::cout << "t4_reference:" << std::endl << t_ref << std::endl; std::cout << "t4_rotator:" << std::endl << t_rotator << std::endl; } Vec_t v_back{Fix::rotator.rotate_back(v_rotator)}; Mat_t m_back{Fix::rotator.rotate_back(m_rotator)}; Ten_t t_back{Fix::rotator.rotate_back(t_rotator)}; v_error = (v_back - v).norm() / v.norm(); BOOST_CHECK_LT(v_error, tol); m_error = (m_back - m).norm() / m.norm(); BOOST_CHECK_LT(m_error, tol); t_error = (t_back - t).norm() / t.norm(); BOOST_CHECK_LT(t_error, tol); } /* ---------------------------------------------------------------------- */ using threeD_list = boost::mpl::list>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(rotation_matrix_test, Fix, threeD_list, Fix) { using Mat_t = typename Fix::Mat_t; auto c{Eigen::cos(Fix::euler.array())}; Real c_1{c[0]}, c_2{c[1]}, c_3{c[2]}; auto s{Eigen::sin(Fix::euler.array())}; Real s_1{s[0]}, s_2{s[1]}, s_3{s[2]}; Mat_t rot_ref; switch (Fix::EulerOrder) { case RotationOrder::ZXYTaitBryan: { rot_ref << c_1 * c_3 - s_1 * s_2 * s_3, -c_2 * s_1, c_1 * s_3 + c_3 * s_1 * s_2, c_3 * s_1 + c_1 * s_2 * s_3, c_1 * c_2, s_1 * s_3 - c_1 * c_3 * s_2, -c_2 * s_3, s_2, c_2 * c_3; break; } default: { BOOST_CHECK(false); break; } } auto err{(rot_ref - Fix::rotator.get_rot_mat()).norm()}; BOOST_CHECK_LT(err, tol); if (err >= tol) { std::cout << "Reference:" << std::endl << rot_ref << std::endl; std::cout << "Rotator:" << std::endl << Fix::rotator.get_rot_mat() << std::endl; } } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_goodies.hh b/tests/test_goodies.hh index 4bc39aa..3afdcf3 100644 --- a/tests/test_goodies.hh +++ b/tests/test_goodies.hh @@ -1,229 +1,228 @@ /** * @file test_goodies.hh * * @author Till Junge * * @date 27 Sep 2017 * * @brief helpers for testing * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef TESTS_TEST_GOODIES_HH_ #define TESTS_TEST_GOODIES_HH_ #include #include #include #include #include namespace muGrid { namespace testGoodies { - + template struct dimFixture { constexpr static Dim_t dim{Dim}; }; using dimlist = boost::mpl::list, dimFixture, dimFixture>; /* ---------------------------------------------------------------------- */ template class RandRange { public: RandRange() : rd(), gen(rd()) {} template std::enable_if_t::value, dummyT> randval(T && lower, T && upper) { static_assert(std::is_same::value, "SFINAE"); auto distro = std::uniform_real_distribution(lower, upper); return distro(this->gen); } template std::enable_if_t::value, dummyT> randval(T && lower, T && upper) { static_assert(std::is_same::value, "SFINAE"); auto distro = std::uniform_int_distribution(lower, upper); return distro(this->gen); } private: std::random_device rd; std::default_random_engine gen; }; /** * explicit computation of linearisation of PK1 stress for an * objective Hooke's law. This implementation is not meant to be * efficient, but te reflect exactly the formulation in Curnier * 2000, "Méthodes numériques en mécanique des solides" for * reference and testing */ template decltype(auto) objective_hooke_explicit(Real lambda, Real mu, const Matrices::Tens2_t & F) { using Matrices::Tens2_t; using Matrices::Tens4_t; using Matrices::tensmult; using T2 = Tens2_t; using T4 = Tens4_t; T2 P; T2 I = P.Identity(); T4 K; // See Curnier, 2000, "Méthodes numériques en mécanique des // solides", p 252, (6.95b) Real Fjrjr = (F.array() * F.array()).sum(); T2 Fjrjm = F.transpose() * F; P.setZero(); for (Dim_t i = 0; i < Dim; ++i) { for (Dim_t m = 0; m < Dim; ++m) { P(i, m) += lambda / 2 * (Fjrjr - Dim) * F(i, m); for (Dim_t r = 0; r < Dim; ++r) { P(i, m) += mu * F(i, r) * (Fjrjm(r, m) - I(r, m)); } } } // See Curnier, 2000, "Méthodes numériques en mécanique des solides", p // 252 Real Fkrkr = (F.array() * F.array()).sum(); T2 Fkmkn = F.transpose() * F; T2 Fisjs = F * F.transpose(); K.setZero(); for (Dim_t i = 0; i < Dim; ++i) { for (Dim_t j = 0; j < Dim; ++j) { for (Dim_t m = 0; m < Dim; ++m) { for (Dim_t n = 0; n < Dim; ++n) { get(K, i, m, j, n) = (lambda * ((Fkrkr - Dim) / 2 * I(i, j) * I(m, n) + F(i, m) * F(j, n)) + mu * (I(i, j) * Fkmkn(m, n) + Fisjs(i, j) * I(m, n) - I(i, j) * I(m, n) + F(i, n) * F(j, m))); } } } } return std::make_tuple(P, K); } /** * takes a 4th-rank tensor and returns a copy with the last two * dimensions switched. This is particularly useful to check for * identity between a stiffness tensors computed the regular way * and the Geers way. For testing, not efficient. */ template inline auto right_transpose(const Eigen::MatrixBase & t4) -> Eigen::Matrix { - using muGrid::get; constexpr Dim_t Rows{Derived::RowsAtCompileTime}; constexpr Dim_t Cols{Derived::ColsAtCompileTime}; constexpr Dim_t DimSq{Rows}; using T = typename Derived::Scalar; static_assert(Rows == Cols, "only square problems"); - constexpr Dim_t Dim{muGrid::ct_sqrt(DimSq)}; + constexpr Dim_t Dim{ct_sqrt(DimSq)}; static_assert((Dim == twoD) or (Dim == threeD), "only two- and three-dimensional problems"); - static_assert(muGrid::ipow(Dim, 2) == DimSq, + static_assert(ipow(Dim, 2) == DimSq, "The array is not a valid fourth order tensor"); - auto retval{muGrid::T4Mat::Zero()}; + T4Mat retval{T4Mat::Zero()}; /** * Note: this looks like it's doing a left transpose, but in * reality, np is rowmajor in all directions, so from numpy to * eigen, we need to transpose left, center, and * right. Therefore, if we want to right-transpose, then we need * to transpose once left, once center, and *twice* right, which * is equal to just left and center transpose. Center-transpose * happens automatically when eigen parses the input (which is * naturally written in rowmajor but interpreted into colmajor), * so all that's left to do is (ironically) the subsequent * left-transpose */ for (int i{0}; i < Dim; ++i) { for (int j{0}; j < Dim; ++j) { for (int k{0}; k < Dim; ++k) { for (int l{0}; l < Dim; ++l) { get(retval, i, j, k, l) = get(t4, j, i, k, l); } } } } return retval; } /** * recomputes a colmajor representation of a fourth-rank rowmajor * tensor. this is useful when comparing to reference results * computed in numpy */ template inline auto from_numpy(const Eigen::MatrixBase & t4_np) -> Eigen::Matrix { constexpr Dim_t Rows{Derived::RowsAtCompileTime}; constexpr Dim_t Cols{Derived::ColsAtCompileTime}; constexpr Dim_t DimSq{Rows}; using T = typename Derived::Scalar; static_assert(Rows == Cols, "only square problems"); - constexpr Dim_t Dim{muGrid::ct_sqrt(DimSq)}; + constexpr Dim_t Dim{ct_sqrt(DimSq)}; static_assert((Dim == twoD) or (Dim == threeD), "only two- and three-dimensional problems"); - static_assert(muGrid::ipow(Dim, 2) == DimSq, + static_assert(ipow(Dim, 2) == DimSq, "The array is not a valid fourth order tensor"); - auto retval{muGrid::T4Mat::Zero()}; - auto intermediate{muGrid::T4Mat::Zero()}; + T4Mat retval{T4Mat::Zero()}; + T4Mat intermediate{T4Mat::Zero()}; // transpose rows for (int row{0}; row < DimSq; ++row) { for (int i{0}; i < Dim; ++i) { for (int j{0}; j < Dim; ++j) { intermediate(row, i + Dim * j) = t4_np(row, i * Dim + j); } } } // transpose columns for (int col{0}; col < DimSq; ++col) { for (int i{0}; i < Dim; ++i) { for (int j{0}; j < Dim; ++j) { retval(i + Dim * j, col) = intermediate(i * Dim + j, col); } } } return retval; } } // namespace testGoodies } // namespace muGrid #endif // TESTS_TEST_GOODIES_HH_ diff --git a/tests/test_material_evaluator.cc b/tests/test_material_evaluator.cc index 44630ea..3429adb 100644 --- a/tests/test_material_evaluator.cc +++ b/tests/test_material_evaluator.cc @@ -1,288 +1,287 @@ /** * @file test_material_evaluator.cc * * @author Till Junge * * @date 13 Jan 2019 * * @brief tests for the material evaluator mechanism * * Copyright © 2019 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "tests.hh" -#include "common/T4_map_proxy.hh" #include "materials/material_linear_elastic2.hh" #include "materials/material_evaluator.hh" +#include + #include "Eigen/Dense" namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_evaluator_tests); /* ---------------------------------------------------------------------- */ BOOST_AUTO_TEST_CASE(without_per_pixel_data) { using Mat_t = MaterialLinearElastic1; constexpr Real Young{210e9}; constexpr Real Poisson{.33}; auto mat_eval = Mat_t::make_evaluator(Young, Poisson); auto & mat = *std::get<0>(mat_eval); auto & evaluator = std::get<1>(mat_eval); using T2_t = Eigen::Matrix; - using T4_t = T4Mat; + using T4_t = muGrid::T4Mat; const T2_t F{(T2_t::Random() - (T2_t::Ones() * .5)) * 1e-4 + T2_t::Identity()}; const T2_t eps{ .5 * ((F - T2_t::Identity()) + (F - T2_t::Identity()).transpose())}; /* * at this point, the evaluator has been created, but the underlying * material still has zero pixels. Evaluation is not yet possible, and * trying to do so has to fail with an explicit error message */ BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain), std::runtime_error); mat.add_pixel({}); const T2_t sigma{evaluator.evaluate_stress(eps, Formulation::small_strain)}; const T2_t P{evaluator.evaluate_stress(F, Formulation::finite_strain)}; auto J{F.determinant()}; auto P_reconstruct{J * sigma * F.inverse().transpose()}; auto error_comp{[](const auto & a, const auto & b) { return (a - b).norm() / (a + b).norm(); }}; auto error{error_comp(P, P_reconstruct)}; constexpr Real small_strain_tol{1e-3}; if (not(error <= small_strain_tol)) { std::cout << "F =" << std::endl << F << std::endl; std::cout << "ε =" << std::endl << eps << std::endl; std::cout << "P =" << std::endl << P << std::endl; std::cout << "σ =" << std::endl << sigma << std::endl; std::cout << "P_reconstructed =" << std::endl << P_reconstruct << std::endl; } BOOST_CHECK_LE(error, small_strain_tol); T2_t sigma2, P2; T4_t C, K; std::tie(sigma2, C) = - evaluator.evaluate_stress_tangent(eps, Formulation::small_strain); + evaluator.evaluate_stress_tangent(eps, Formulation::small_strain); std::tie(P2, K) = - evaluator.evaluate_stress_tangent(F, Formulation::finite_strain); + evaluator.evaluate_stress_tangent(F, Formulation::finite_strain); error = error_comp(sigma2, sigma); BOOST_CHECK_LE(error, tol); error = error_comp(P2, P); BOOST_CHECK_LE(error, tol); error = error_comp(C, K); if (not(error <= small_strain_tol)) { std::cout << "F =" << std::endl << F << std::endl; std::cout << "ε =" << std::endl << eps << std::endl; std::cout << "P =" << std::endl << P << std::endl; std::cout << "σ =" << std::endl << sigma << std::endl; std::cout << "K =" << std::endl << K << std::endl; std::cout << "C =" << std::endl << C << std::endl; } BOOST_CHECK_LE(error, small_strain_tol); - mat.add_pixel({1}); /* * Now, the material has two pixels, and evaluating it would be ambiguous. * It should fail with an explicit error message */ BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain), std::runtime_error); } /* ---------------------------------------------------------------------- */ BOOST_AUTO_TEST_CASE(with_per_pixel_data) { using Mat_t = MaterialLinearElastic2; constexpr Real Young{210e9}; constexpr Real Poisson{.33}; auto mat_eval{Mat_t::make_evaluator(Young, Poisson)}; auto & mat{*std::get<0>(mat_eval)}; auto & evaluator{std::get<1>(mat_eval)}; - using T2_t = Eigen::Matrix; - using T4_t = T4Mat; + using T4_t = muGrid::T4Mat; const T2_t F{(T2_t::Random() - (T2_t::Ones() * .5)) * 1e-4 + T2_t::Identity()}; const T2_t eps{ .5 * ((F - T2_t::Identity()) + (F - T2_t::Identity()).transpose())}; BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain), std::runtime_error); T2_t eigen_strain{[](auto x) { return 1e-4 * (x + x.transpose()); }(T2_t::Random() - T2_t::Ones() * .5)}; mat.add_pixel({}, eigen_strain); const T2_t sigma{evaluator.evaluate_stress(eps, Formulation::small_strain)}; const T2_t P{evaluator.evaluate_stress(F, Formulation::finite_strain)}; auto J{F.determinant()}; auto P_reconstruct{J * sigma * F.inverse().transpose()}; auto error_comp{[](const auto & a, const auto & b) { return (a - b).norm() / (a + b).norm(); }}; auto error{error_comp(P, P_reconstruct)}; constexpr Real small_strain_tol{1e-3}; if (not(error <= small_strain_tol)) { std::cout << "F =" << std::endl << F << std::endl; std::cout << "ε =" << std::endl << eps << std::endl; std::cout << "P =" << std::endl << P << std::endl; std::cout << "σ =" << std::endl << sigma << std::endl; std::cout << "P_reconstructed =" << std::endl << P_reconstruct << std::endl; } BOOST_CHECK_LE(error, small_strain_tol); T2_t sigma2, P2; T4_t C, K; std::tie(sigma2, C) = - evaluator.evaluate_stress_tangent(eps, Formulation::small_strain); + evaluator.evaluate_stress_tangent(eps, Formulation::small_strain); std::tie(P2, K) = - evaluator.evaluate_stress_tangent(F, Formulation::finite_strain); + evaluator.evaluate_stress_tangent(F, Formulation::finite_strain); error = error_comp(sigma2, sigma); BOOST_CHECK_LE(error, tol); error = error_comp(P2, P); BOOST_CHECK_LE(error, tol); error = error_comp(C, K); if (not(error <= small_strain_tol)) { std::cout << "F =" << std::endl << F << std::endl; std::cout << "ε =" << std::endl << eps << std::endl; std::cout << "P =" << std::endl << P << std::endl; std::cout << "σ =" << std::endl << sigma << std::endl; std::cout << "K =" << std::endl << K << std::endl; std::cout << "C =" << std::endl << C << std::endl; } BOOST_CHECK_LE(error, small_strain_tol); } /* ---------------------------------------------------------------------- */ BOOST_AUTO_TEST_CASE(tangent_estimation) { using Mat_t = MaterialLinearElastic1; constexpr Real Young{210e9}; constexpr Real Poisson{.33}; auto mat_eval = Mat_t::make_evaluator(Young, Poisson); auto & mat = *std::get<0>(mat_eval); auto & evaluator = std::get<1>(mat_eval); using T2_t = Eigen::Matrix; - using T4_t = T4Mat; + using T4_t = muGrid::T4Mat; const T2_t F{(T2_t::Random() - (T2_t::Ones() * .5)) * 1e-4 + T2_t::Identity()}; const T2_t eps{ .5 * ((F - T2_t::Identity()) + (F - T2_t::Identity()).transpose())}; BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain), std::runtime_error); mat.add_pixel({}); T2_t sigma, P; T4_t C, K; std::tie(sigma, C) = - evaluator.evaluate_stress_tangent(eps, Formulation::small_strain); + evaluator.evaluate_stress_tangent(eps, Formulation::small_strain); std::tie(P, K) = - evaluator.evaluate_stress_tangent(F, Formulation::finite_strain); + evaluator.evaluate_stress_tangent(F, Formulation::finite_strain); constexpr Real linear_step{1.}; constexpr Real nonlin_step{1.e-6}; T4_t C_estim{evaluator.estimate_tangent(eps, Formulation::small_strain, linear_step)}; - T4_t K_estim{evaluator.estimate_tangent(F, Formulation::finite_strain, - nonlin_step)}; + T4_t K_estim{ + evaluator.estimate_tangent(F, Formulation::finite_strain, nonlin_step)}; auto error_comp{[](const auto & a, const auto & b) { return (a - b).norm() / (a + b).norm(); }}; constexpr Real finite_diff_tol{1e-9}; - Real error {error_comp(K, K_estim)}; + Real error{error_comp(K, K_estim)}; if (not(error <= finite_diff_tol)) { std::cout << "K =" << std::endl << K << std::endl; std::cout << "K_estim =" << std::endl << K_estim << std::endl; } BOOST_CHECK_LE(error, finite_diff_tol); error = error_comp(C, C_estim); if (not(error <= tol)) { std::cout << "centred difference:" << std::endl; std::cout << "C =" << std::endl << C << std::endl; std::cout << "C_estim =" << std::endl << C_estim << std::endl; } BOOST_CHECK_LE(error, tol); C_estim = evaluator.estimate_tangent(eps, Formulation::small_strain, linear_step, FiniteDiff::forward); error = error_comp(C, C_estim); if (not(error <= tol)) { std::cout << "forward difference:" << std::endl; std::cout << "C =" << std::endl << C << std::endl; std::cout << "C_estim =" << std::endl << C_estim << std::endl; } BOOST_CHECK_LE(error, tol); C_estim = evaluator.estimate_tangent(eps, Formulation::small_strain, linear_step, FiniteDiff::backward); error = error_comp(C, C_estim); if (not(error <= tol)) { std::cout << "backward difference:" << std::endl; std::cout << "C =" << std::endl << C << std::endl; std::cout << "C_estim =" << std::endl << C_estim << std::endl; } BOOST_CHECK_LE(error, tol); -} + } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_material_hyper_elasto_plastic1.cc b/tests/test_material_hyper_elasto_plastic1.cc index 7c7b281..d1310b7 100644 --- a/tests/test_material_hyper_elasto_plastic1.cc +++ b/tests/test_material_hyper_elasto_plastic1.cc @@ -1,425 +1,429 @@ /** * @file test_material_hyper_elasto_plastic1.cc * * @author Till Junge * * @date 25 Feb 2018 * * @brief Tests for the large-strain Simo-type plastic law implemented * using MaterialMuSpectre * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "boost/mpl/list.hpp" #include "materials/stress_transformations_Kirchhoff.hh" #include "materials/material_hyper_elasto_plastic1.hh" #include "materials/materials_toolbox.hh" #include "tests.hh" #include "test_goodies.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_hyper_elasto_plastic_1); template struct MaterialFixture { using Mat = Mat_t; constexpr static Real K{.833}; // bulk modulus constexpr static Real mu{.386}; // shear modulus constexpr static Real H{.004}; // hardening modulus constexpr static Real tau_y0{.003}; // initial yield stress constexpr static Real young{MatTB::convert_elastic_modulus< ElasticModulus::Young, ElasticModulus::Bulk, ElasticModulus::Shear>( K, mu)}; constexpr static Real poisson{MatTB::convert_elastic_modulus< ElasticModulus::Poisson, ElasticModulus::Bulk, ElasticModulus::Shear>( K, mu)}; MaterialFixture() : mat("Name", young, poisson, tau_y0, H) {} constexpr static Dim_t sdim{Mat_t::sdim()}; constexpr static Dim_t mdim{Mat_t::mdim()}; Mat_t mat; }; using mats = boost::mpl::list< MaterialFixture>, MaterialFixture>, MaterialFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constructor, Fix, mats, Fix) { BOOST_CHECK_EQUAL("Name", Fix::mat.get_name()); auto & mat{Fix::mat}; auto sdim{Fix::sdim}; auto mdim{Fix::mdim}; BOOST_CHECK_EQUAL(sdim, mat.sdim()); BOOST_CHECK_EQUAL(mdim, mat.mdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_stress, Fix, mats, Fix) { // This test uses precomputed reference values (computed using // elasto-plasticity.py) for the 3d case only // need higher tol because of printout precision of reference solutions constexpr Real hi_tol{1e-8}; constexpr Dim_t mdim{Fix::mdim}, sdim{Fix::sdim}; constexpr bool has_precomputed_values{(mdim == sdim) && (mdim == threeD)}; constexpr bool verbose{false}; using Strain_t = Eigen::Matrix; using traits = MaterialMuSpectre_traits>; using LColl_t = typename traits::LFieldColl_t; - using StrainStField_t = - StateField>; - using FlowStField_t = StateField>; + using StrainStField_t = muGrid::StateField< + muGrid::TensorField>; + using FlowStField_t = + muGrid::StateField>; // using StrainStRef_t = typename traits::LStrainMap_t::reference; // using ScalarStRef_t = typename traits::LScalarMap_t::reference; // create statefields LColl_t coll{}; coll.add_pixel({0}); coll.initialise(); - auto & F_{make_statefield("previous gradient", coll)}; - auto & be_{ - make_statefield("previous elastic strain", coll)}; - auto & eps_{make_statefield("plastic flow", coll)}; + auto & F_{ + muGrid::make_statefield("previous gradient", coll)}; + auto & be_{muGrid::make_statefield( + "previous elastic strain", coll)}; + auto & eps_{muGrid::make_statefield("plastic flow", coll)}; auto F_prev{F_.get_map()}; F_prev[0].current() = Strain_t::Identity(); auto be_prev{be_.get_map()}; be_prev[0].current() = Strain_t::Identity(); auto eps_prev{eps_.get_map()}; eps_prev[0].current() = 0; // elastic deformation Strain_t F{Strain_t::Identity()}; F(0, 1) = 1e-5; F_.cycle(); be_.cycle(); eps_.cycle(); Strain_t stress{ Fix::mat.evaluate_stress(F, F_prev[0], be_prev[0], eps_prev[0])}; if (has_precomputed_values) { Strain_t tau_ref{}; tau_ref << 1.92999522e-11, 3.86000000e-06, 0.00000000e+00, 3.86000000e-06, -1.93000510e-11, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -2.95741950e-17; Real error{(tau_ref - stress).norm() / tau_ref.norm()}; BOOST_CHECK_LT(error, hi_tol); Strain_t be_ref{}; be_ref << 1.00000000e+00, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00; error = (be_ref - be_prev[0].current()).norm() / be_ref.norm(); BOOST_CHECK_LT(error, hi_tol); Real ep_ref{0}; error = ep_ref - eps_prev[0].current(); BOOST_CHECK_LT(error, hi_tol); } if (verbose) { std::cout << "τ =" << std::endl << stress << std::endl; std::cout << "F =" << std::endl << F << std::endl; std::cout << "Fₜ =" << std::endl << F_prev[0].current() << std::endl; std::cout << "bₑ =" << std::endl << be_prev[0].current() << std::endl; std::cout << "εₚ =" << std::endl << eps_prev[0].current() << std::endl; } F_.cycle(); be_.cycle(); eps_.cycle(); // plastic deformation F(0, 1) = .2; stress = Fix::mat.evaluate_stress(F, F_prev[0], be_prev[0], eps_prev[0]); if (has_precomputed_values) { Strain_t tau_ref{}; tau_ref << 1.98151335e-04, 1.98151335e-03, 0.00000000e+00, 1.98151335e-03, -1.98151335e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.60615155e-16; Real error{(tau_ref - stress).norm() / tau_ref.norm()}; BOOST_CHECK_LT(error, hi_tol); Strain_t be_ref{}; be_ref << 1.00052666, 0.00513348, 0., 0.00513348, 0.99949996, 0., 0., 0., 1.; error = (be_ref - be_prev[0].current()).norm() / be_ref.norm(); BOOST_CHECK_LT(error, hi_tol); Real ep_ref{0.11229988}; error = (ep_ref - eps_prev[0].current()) / ep_ref; BOOST_CHECK_LT(error, hi_tol); } if (verbose) { std::cout << "Post Cycle" << std::endl; std::cout << "τ =" << std::endl << stress << std::endl << "F =" << std::endl << F << std::endl << "Fₜ =" << std::endl << F_prev[0].current() << std::endl << "bₑ =" << std::endl << be_prev[0].current() << std::endl << "εₚ =" << std::endl << eps_prev[0].current() << std::endl; } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_stiffness, Fix, mats, Fix) { // This test uses precomputed reference values (computed using // elasto-plasticity.py) for the 3d case only // need higher tol because of printout precision of reference solutions constexpr Real hi_tol{2e-7}; constexpr Dim_t mdim{Fix::mdim}, sdim{Fix::sdim}; constexpr bool has_precomputed_values{(mdim == sdim) && (mdim == threeD)}; constexpr bool verbose{has_precomputed_values && false}; using Strain_t = Eigen::Matrix; - using Stiffness_t = T4Mat; + using Stiffness_t = muGrid::T4Mat; using traits = MaterialMuSpectre_traits>; using LColl_t = typename traits::LFieldColl_t; - using StrainStField_t = - StateField>; - using FlowStField_t = StateField>; + using StrainStField_t = muGrid::StateField< + muGrid::TensorField>; + using FlowStField_t = + muGrid::StateField>; // using StrainStRef_t = typename traits::LStrainMap_t::reference; // using ScalarStRef_t = typename traits::LScalarMap_t::reference; // create statefields LColl_t coll{}; coll.add_pixel({0}); coll.initialise(); - auto & F_{make_statefield("previous gradient", coll)}; - auto & be_{ - make_statefield("previous elastic strain", coll)}; - auto & eps_{make_statefield("plastic flow", coll)}; + auto & F_{ + muGrid::make_statefield("previous gradient", coll)}; + auto & be_{muGrid::make_statefield( + "previous elastic strain", coll)}; + auto & eps_{muGrid::make_statefield("plastic flow", coll)}; auto F_prev{F_.get_map()}; F_prev[0].current() = Strain_t::Identity(); auto be_prev{be_.get_map()}; be_prev[0].current() = Strain_t::Identity(); auto eps_prev{eps_.get_map()}; eps_prev[0].current() = 0; // elastic deformation Strain_t F{Strain_t::Identity()}; F(0, 1) = 1e-5; F_.cycle(); be_.cycle(); eps_.cycle(); Strain_t stress{}; Stiffness_t stiffness{}; std::tie(stress, stiffness) = Fix::mat.evaluate_stress_tangent(F, F_prev[0], be_prev[0], eps_prev[0]); if (has_precomputed_values) { Strain_t tau_ref{}; tau_ref << 1.92999522e-11, 3.86000000e-06, 0.00000000e+00, 3.86000000e-06, -1.93000510e-11, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -2.95741950e-17; Real error{(tau_ref - stress).norm() / tau_ref.norm()}; BOOST_CHECK_LT(error, hi_tol); Strain_t be_ref{}; be_ref << 1.00000000e+00, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00; error = (be_ref - be_prev[0].current()).norm() / be_ref.norm(); BOOST_CHECK_LT(error, hi_tol); Real ep_ref{0}; error = ep_ref - eps_prev[0].current(); BOOST_CHECK_LT(error, hi_tol); Stiffness_t temp; temp << 1.34766667e+00, 3.86000000e-06, 0.00000000e+00, -3.86000000e-06, 5.75666667e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 5.75666667e-01, -3.61540123e-17, 3.86000000e-01, 0.00000000e+00, 3.86000000e-01, 7.12911684e-17, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 3.86000000e-01, 0.00000000e+00, 0.00000000e+00, -1.93000000e-06, 3.86000000e-01, 1.93000000e-06, 0.00000000e+00, -3.61540123e-17, 3.86000000e-01, 0.00000000e+00, 3.86000000e-01, 7.12911684e-17, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 5.75666667e-01, -3.86000000e-06, 0.00000000e+00, 3.86000000e-06, 1.34766667e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 5.75666667e-01, 0.00000000e+00, 0.00000000e+00, -1.93000000e-06, 0.00000000e+00, 0.00000000e+00, 3.86000000e-01, 1.93000000e-06, 3.86000000e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 3.86000000e-01, 0.00000000e+00, 0.00000000e+00, -1.93000000e-06, 3.86000000e-01, 1.93000000e-06, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -1.93000000e-06, 0.00000000e+00, 0.00000000e+00, 3.86000000e-01, 1.93000000e-06, 3.86000000e-01, 0.00000000e+00, 5.75666667e-01, 2.61999996e-17, 0.00000000e+00, 2.61999996e-17, 5.75666667e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.34766667e+00; - Stiffness_t K4b_ref{testGoodies::from_numpy(temp)}; + Stiffness_t K4b_ref{muGrid::testGoodies::from_numpy(temp)}; error = (K4b_ref - stiffness).norm() / K4b_ref.norm(); BOOST_CHECK_LT(error, hi_tol); if (not(error < hi_tol)) { std::cout << "stiffness reference:\n" << K4b_ref << std::endl; std::cout << "stiffness computed:\n" << stiffness << std::endl; } } if (verbose) { std::cout << "C₄ =" << std::endl << stiffness << std::endl; } F_.cycle(); be_.cycle(); eps_.cycle(); // plastic deformation F(0, 1) = .2; std::tie(stress, stiffness) = Fix::mat.evaluate_stress_tangent(F, F_prev[0], be_prev[0], eps_prev[0]); if (has_precomputed_values) { Strain_t tau_ref{}; tau_ref << 1.98151335e-04, 1.98151335e-03, 0.00000000e+00, 1.98151335e-03, -1.98151335e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.60615155e-16; Real error{(tau_ref - stress).norm() / tau_ref.norm()}; BOOST_CHECK_LT(error, hi_tol); Strain_t be_ref{}; be_ref << 1.00052666, 0.00513348, 0., 0.00513348, 0.99949996, 0., 0., 0., 1.; error = (be_ref - be_prev[0].current()).norm() / be_ref.norm(); BOOST_CHECK_LT(error, hi_tol); Real ep_ref{0.11229988}; error = (ep_ref - eps_prev[0].current()) / ep_ref; BOOST_CHECK_LT(error, hi_tol); Stiffness_t temp{}; temp << 8.46343327e-01, 1.11250597e-03, 0.00000000e+00, -2.85052074e-03, 8.26305692e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 8.26350980e-01, -8.69007382e-04, 1.21749295e-03, 0.00000000e+00, 1.61379562e-03, 8.69007382e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 5.58242059e-18, 0.00000000e+00, 0.00000000e+00, 9.90756677e-03, 0.00000000e+00, 0.00000000e+00, -9.90756677e-04, 1.01057181e-02, 9.90756677e-04, 0.00000000e+00, -8.69007382e-04, 1.21749295e-03, 0.00000000e+00, 1.61379562e-03, 8.69007382e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 5.58242059e-18, 8.26305692e-01, -1.11250597e-03, 0.00000000e+00, 2.85052074e-03, 8.46343327e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 8.26350980e-01, 0.00000000e+00, 0.00000000e+00, -9.90756677e-04, 0.00000000e+00, 0.00000000e+00, 1.01057181e-02, 9.90756677e-04, 9.90756677e-03, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 9.90756677e-03, 0.00000000e+00, 0.00000000e+00, -9.90756677e-04, 1.01057181e-02, 9.90756677e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -9.90756677e-04, 0.00000000e+00, 0.00000000e+00, 1.01057181e-02, 9.90756677e-04, 9.90756677e-03, 0.00000000e+00, 8.26350980e-01, 0.00000000e+00, 0.00000000e+00, 1.38777878e-17, 8.26350980e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 8.46298039e-01; - Stiffness_t K4b_ref{testGoodies::from_numpy(temp)}; + Stiffness_t K4b_ref{muGrid::testGoodies::from_numpy(temp)}; error = (K4b_ref - stiffness).norm() / K4b_ref.norm(); error = (K4b_ref - stiffness).norm() / K4b_ref.norm(); BOOST_CHECK_LT(error, hi_tol); if (not(error < hi_tol)) { std::cout << "stiffness reference:\n" << K4b_ref << std::endl; std::cout << "stiffness computed:\n" << stiffness << std::endl; } // check also whether pull_back is correct Stiffness_t intermediate{stiffness}; Stiffness_t zero_mediate{Stiffness_t::Zero()}; for (int i{0}; i < mdim; ++i) { for (int j{0}; j < mdim; ++j) { for (int m{0}; m < mdim; ++m) { const auto & k{i}; const auto & l{j}; // k,m inverted for right transpose - get(zero_mediate, i, j, k, m) -= stress(l, m); - get(intermediate, i, j, k, m) -= stress(l, m); + muGrid::get(zero_mediate, i, j, k, m) -= stress(l, m); + muGrid::get(intermediate, i, j, k, m) -= stress(l, m); } } } temp << 8.46145176e-01, -8.69007382e-04, 0.00000000e+00, -2.85052074e-03, 8.26305692e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 8.26350980e-01, -2.85052074e-03, 1.41564428e-03, 0.00000000e+00, 1.61379562e-03, 8.69007382e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 5.58242059e-18, 0.00000000e+00, 0.00000000e+00, 9.90756677e-03, 0.00000000e+00, 0.00000000e+00, -9.90756677e-04, 1.01057181e-02, 9.90756677e-04, 0.00000000e+00, -8.69007382e-04, 1.21749295e-03, 0.00000000e+00, 1.41564428e-03, -1.11250597e-03, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 5.58242059e-18, 8.26305692e-01, -1.11250597e-03, 0.00000000e+00, 8.69007382e-04, 8.46541479e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 8.26350980e-01, 0.00000000e+00, 0.00000000e+00, -9.90756677e-04, 0.00000000e+00, 0.00000000e+00, 1.01057181e-02, 9.90756677e-04, 9.90756677e-03, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 9.90756677e-03, 0.00000000e+00, 0.00000000e+00, -9.90756677e-04, 9.90756677e-03, -9.90756677e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -9.90756677e-04, 0.00000000e+00, 0.00000000e+00, 1.01057181e-02, -9.90756677e-04, 1.01057181e-02, 0.00000000e+00, 8.26350980e-01, 0.00000000e+00, 0.00000000e+00, 1.38777878e-17, 8.26350980e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 8.46298039e-01; - Stiffness_t K4c_ref{testGoodies::from_numpy(temp)}; + Stiffness_t K4c_ref{muGrid::testGoodies::from_numpy(temp)}; error = (K4b_ref + zero_mediate - K4c_ref).norm() / zero_mediate.norm(); BOOST_CHECK_LT(error, hi_tol); // rel error on small difference between // inexacly read doubles if (not(error < hi_tol)) { std::cout << "decrement reference:\n" << K4c_ref - K4b_ref << std::endl; std::cout << "zero_mediate computed:\n" << zero_mediate << std::endl; } error = (K4c_ref - intermediate).norm() / K4c_ref.norm(); BOOST_CHECK_LT(error, hi_tol); if (not(error < hi_tol)) { std::cout << "stiffness reference:\n" << K4c_ref << std::endl; std::cout << "stiffness computed:\n" << intermediate << std::endl; std::cout << "zero-mediate computed:\n" << zero_mediate << std::endl; std::cout << "difference:\n" << K4c_ref - intermediate << std::endl; } } if (verbose) { std::cout << "Post Cycle" << std::endl; std::cout << "C₄ =" << std::endl << stiffness << std::endl; } } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(stress_strain_test, Fix, mats, Fix) {} BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_material_linear_elastic1.cc b/tests/test_material_linear_elastic1.cc index ed36e30..c2ad947 100644 --- a/tests/test_material_linear_elastic1.cc +++ b/tests/test_material_linear_elastic1.cc @@ -1,235 +1,236 @@ /** * @file test_material_linear_elastic1.cc * * @author Till Junge * * @date 28 Nov 2017 * * @brief Tests for the large-strain, objective Hooke's law, implemented in * the convenient strategy (i.e., using MaterialMuSpectre), also used * to test parts of MaterialLinearElastic2 * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include #include #include "materials/material_linear_elastic1.hh" #include "materials/material_linear_elastic2.hh" #include "tests.hh" #include "tests/test_goodies.hh" -#include "common/field_collection.hh" -#include "common/iterators.hh" +#include +#include namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_linear_elastic_1); template struct MaterialFixture { using Mat = Mat_t; constexpr static Real lambda{2}, mu{1.5}; constexpr static Real young{mu * (3 * lambda + 2 * mu) / (lambda + mu)}; constexpr static Real poisson{lambda / (2 * (lambda + mu))}; MaterialFixture() : mat("Name", young, poisson) {} constexpr static Dim_t sdim{Mat_t::sdim()}; constexpr static Dim_t mdim{Mat_t::mdim()}; Mat_t mat; }; template struct has_internals { constexpr static bool value{false}; }; template struct has_internals> { constexpr static bool value{true}; }; using mat_list = boost::mpl::list>, MaterialFixture>, MaterialFixture>, MaterialFixture>, MaterialFixture>, MaterialFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constructor, Fix, mat_list, Fix) { BOOST_CHECK_EQUAL("Name", Fix::mat.get_name()); auto & mat{Fix::mat}; auto sdim{Fix::sdim}; auto mdim{Fix::mdim}; BOOST_CHECK_EQUAL(sdim, mat.sdim()); BOOST_CHECK_EQUAL(mdim, mat.mdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_add_pixel, Fix, mat_list, Fix) { auto & mat{Fix::mat}; constexpr Dim_t sdim{Fix::sdim}; - testGoodies::RandRange rng; + muGrid::testGoodies::RandRange rng; const Dim_t nb_pixel{7}, box_size{17}; using Ccoord = Ccoord_t; for (Dim_t i = 0; i < nb_pixel; ++i) { Ccoord c; for (Dim_t j = 0; j < sdim; ++j) { c[j] = rng.randval(0, box_size); } if (!has_internals::value) { BOOST_CHECK_NO_THROW(mat.add_pixel(c)); } } BOOST_CHECK_NO_THROW(mat.initialise()); } template struct MaterialFixtureFilled : public MaterialFixture { using Mat = Mat_t; constexpr static Dim_t box_size{3}; MaterialFixtureFilled() : MaterialFixture() { using Ccoord = Ccoord_t; - Ccoord cube{CcoordOps::get_cube(box_size)}; - CcoordOps::Pixels pixels(cube); + Ccoord cube{muGrid::CcoordOps::get_cube(box_size)}; + muGrid::CcoordOps::Pixels pixels(cube); for (auto pixel : pixels) { this->mat.add_pixel(pixel); } this->mat.initialise(); } }; using mat_fill = boost::mpl::list< MaterialFixtureFilled>, MaterialFixtureFilled>, MaterialFixtureFilled>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_law, Fix, mat_fill, Fix) { - constexpr auto cube{CcoordOps::get_cube(Fix::box_size)}; - constexpr auto loc{CcoordOps::get_cube(0)}; + constexpr auto cube{muGrid::CcoordOps::get_cube(Fix::box_size)}; + constexpr auto loc{muGrid::CcoordOps::get_cube(0)}; auto & mat{Fix::mat}; - using FC_t = GlobalFieldCollection; + using FC_t = muGrid::GlobalFieldCollection; FC_t globalfields; - auto & F{make_field( + auto & F{muGrid::make_field( "Transformation Gradient", globalfields)}; - auto & P1 = make_field( + auto & P1 = muGrid::make_field( "Nominal Stress1", globalfields); // to be computed alone - auto & P2 = make_field( + auto & P2 = muGrid::make_field( "Nominal Stress2", globalfields); // to be computed with tangent - auto & K = make_field( + auto & K = muGrid::make_field( "Tangent Moduli", globalfields); // to be computed with tangent - auto & Pr = make_field( + auto & Pr = muGrid::make_field( "Nominal Stress reference", globalfields); - auto & Kr = make_field( + auto & Kr = muGrid::make_field( "Tangent Moduli reference", globalfields); // to be computed with tangent globalfields.initialise(cube, loc); static_assert( std::is_same::value, "oh oh"); static_assert( std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert( std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); using traits = MaterialMuSpectre_traits; { // block to contain not-constant gradient map typename traits::StressMap_t grad_map( globalfields["Transformation Gradient"]); for (auto F_ : grad_map) { F_.setRandom(); } grad_map[0] = grad_map[0].Identity(); // identifiable gradients for debug grad_map[1] = 1.2 * grad_map[1].Identity(); // ditto } // compute stresses using material mat.compute_stresses(globalfields["Transformation Gradient"], globalfields["Nominal Stress1"], Formulation::finite_strain); // compute stresses and tangent moduli using material BOOST_CHECK_THROW( mat.compute_stresses_tangent(globalfields["Transformation Gradient"], globalfields["Nominal Stress2"], globalfields["Nominal Stress2"], Formulation::finite_strain), std::runtime_error); mat.compute_stresses_tangent(globalfields["Transformation Gradient"], globalfields["Nominal Stress2"], globalfields["Tangent Moduli"], Formulation::finite_strain); typename traits::StrainMap_t Fmap(globalfields["Transformation Gradient"]); typename traits::StressMap_t Pmap_ref( globalfields["Nominal Stress reference"]); typename traits::TangentMap_t Kmap_ref( globalfields["Tangent Moduli reference"]); for (auto tup : akantu::zip(Fmap, Pmap_ref, Kmap_ref)) { auto F_ = std::get<0>(tup); auto P_ = std::get<1>(tup); auto K_ = std::get<2>(tup); - std::tie(P_, K_) = testGoodies::objective_hooke_explicit( - Fix::lambda, Fix::mu, F_); + std::tie(P_, K_) = + muGrid::testGoodies::objective_hooke_explicit(Fix::lambda, + Fix::mu, F_); } typename traits::StressMap_t Pmap_1(globalfields["Nominal Stress1"]); for (auto tup : akantu::zip(Pmap_ref, Pmap_1)) { auto P_r = std::get<0>(tup); auto P_1 = std::get<1>(tup); Real error = (P_r - P_1).norm(); BOOST_CHECK_LT(error, tol); } typename traits::StressMap_t Pmap_2(globalfields["Nominal Stress2"]); typename traits::TangentMap_t Kmap(globalfields["Tangent Moduli"]); for (auto tup : akantu::zip(Pmap_ref, Pmap_2, Kmap_ref, Kmap)) { auto P_r = std::get<0>(tup); auto P = std::get<1>(tup); Real error = (P_r - P).norm(); BOOST_CHECK_LT(error, tol); auto K_r = std::get<2>(tup); auto K = std::get<3>(tup); error = (K_r - K).norm(); BOOST_CHECK_LT(error, tol); } } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_material_linear_elastic2.cc b/tests/test_material_linear_elastic2.cc index 56222ab..420f162 100644 --- a/tests/test_material_linear_elastic2.cc +++ b/tests/test_material_linear_elastic2.cc @@ -1,280 +1,283 @@ /** * @file test_material_linear_elastic2.cc * * @author Till Junge * * @date 04 Feb 2018 * * @brief Tests for the objective Hooke's law with eigenstrains, * (tests that do not require add_pixel are integrated into * `test_material_linear_elastic1.cc` * * @section LICENSE * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include #include #include "materials/material_linear_elastic2.hh" #include "tests.hh" #include "tests/test_goodies.hh" -#include "common/field_collection.hh" -#include "common/iterators.hh" + +#include +#include + namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_linear_elastic_2); template struct MaterialFixture { using Mat = Mat_t; constexpr static Real lambda{2}, mu{1.5}; constexpr static Real young{mu * (3 * lambda + 2 * mu) / (lambda + mu)}; constexpr static Real poisson{lambda / (2 * (lambda + mu))}; MaterialFixture() : mat("Name", young, poisson) {} constexpr static Dim_t sdim{Mat_t::sdim()}; constexpr static Dim_t mdim{Mat_t::mdim()}; Mat_t mat; }; using mat_list = boost::mpl::list>, MaterialFixture>, MaterialFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constructor, Fix, mat_list, Fix) { BOOST_CHECK_EQUAL("Name", Fix::mat.get_name()); auto & mat{Fix::mat}; auto sdim{Fix::sdim}; auto mdim{Fix::mdim}; BOOST_CHECK_EQUAL(sdim, mat.sdim()); BOOST_CHECK_EQUAL(mdim, mat.mdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_add_pixel, Fix, mat_list, Fix) { auto & mat{Fix::mat}; constexpr Dim_t sdim{Fix::sdim}; - testGoodies::RandRange rng; + muGrid::testGoodies::RandRange rng; const Dim_t nb_pixel{7}, box_size{17}; using Ccoord = Ccoord_t; for (Dim_t i = 0; i < nb_pixel; ++i) { Ccoord c; for (Dim_t j = 0; j < sdim; ++j) { c[j] = rng.randval(0, box_size); } Eigen::Matrix Zero = Eigen::Matrix::Zero(); BOOST_CHECK_NO_THROW(mat.add_pixel(c, Zero)); } BOOST_CHECK_NO_THROW(mat.initialise()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_eigenstrain_equivalence, Fix, mat_list, Fix) { auto & mat{Fix::mat}; const Dim_t nb_pixel{2}; - constexpr auto cube{CcoordOps::get_cube(nb_pixel)}; - constexpr auto loc{CcoordOps::get_cube(0)}; + constexpr auto cube{muGrid::CcoordOps::get_cube(nb_pixel)}; + constexpr auto loc{muGrid::CcoordOps::get_cube(0)}; using Mat_t = Eigen::Matrix; - using FC_t = GlobalFieldCollection; + using FC_t = muGrid::GlobalFieldCollection; FC_t globalfields; - auto & F_f{make_field( + auto & F_f{muGrid::make_field( "Transformation Gradient", globalfields)}; - auto & P1_f = make_field( + auto & P1_f = muGrid::make_field( "Nominal Stress1", globalfields); // to be computed alone - auto & K_f = make_field( + auto & K_f = muGrid::make_field( "Tangent Moduli", globalfields); // to be computed with tangent globalfields.initialise(cube, loc); Mat_t zero{Mat_t::Zero()}; Mat_t F{Mat_t::Random() / 100 + Mat_t::Identity()}; Mat_t strain{-.5 * (F + F.transpose()) - Mat_t::Identity()}; using Ccoord = Ccoord_t; Ccoord pix0{0}; Ccoord pix1{1}; mat.add_pixel(pix0, zero); mat.add_pixel(pix1, strain); mat.initialise(); F_f.get_map()[pix0] = -strain; F_f.get_map()[pix1] = zero; mat.compute_stresses_tangent(F_f, P1_f, K_f, Formulation::small_strain); Real error{(P1_f.get_map()[pix0] - P1_f.get_map()[pix1]).norm()}; Real tol{1e-12}; if (error >= tol) { std::cout << "error = " << error << " >= " << tol << " = tol" << std::endl; std::cout << "P(0) =" << std::endl << P1_f.get_map()[pix0] << std::endl; std::cout << "P(1) =" << std::endl << P1_f.get_map()[pix1] << std::endl; } BOOST_CHECK_LT(error, tol); } template struct MaterialFixtureFilled : public MaterialFixture { using Par = MaterialFixture; using Mat = Mat_t; constexpr static Dim_t box_size{3}; MaterialFixtureFilled() : MaterialFixture() { using Ccoord = Ccoord_t; - Ccoord cube{CcoordOps::get_cube(box_size)}; - CcoordOps::Pixels pixels(cube); + Ccoord cube{muGrid::CcoordOps::get_cube(box_size)}; + muGrid::CcoordOps::Pixels pixels(cube); for (auto pixel : pixels) { Eigen::Matrix Zero = Eigen::Matrix::Zero(); this->mat.add_pixel(pixel, Zero); } this->mat.initialise(); } }; using mat_fill = boost::mpl::list< MaterialFixtureFilled>, MaterialFixtureFilled>, MaterialFixtureFilled>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_law, Fix, mat_fill, Fix) { - constexpr auto cube{CcoordOps::get_cube(Fix::box_size)}; - constexpr auto loc{CcoordOps::get_cube(0)}; + constexpr auto cube{muGrid::CcoordOps::get_cube(Fix::box_size)}; + constexpr auto loc{muGrid::CcoordOps::get_cube(0)}; auto & mat{Fix::mat}; - using FC_t = GlobalFieldCollection; + using FC_t = muGrid::GlobalFieldCollection; FC_t globalfields; - auto & F{make_field( + auto & F{muGrid::make_field( "Transformation Gradient", globalfields)}; - auto & P1 = make_field( + auto & P1 = muGrid::make_field( "Nominal Stress1", globalfields); // to be computed alone - auto & P2 = make_field( + auto & P2 = muGrid::make_field( "Nominal Stress2", globalfields); // to be computed with tangent - auto & K = make_field( + auto & K = muGrid::make_field( "Tangent Moduli", globalfields); // to be computed with tangent - auto & Pr = make_field( + auto & Pr = muGrid::make_field( "Nominal Stress reference", globalfields); - auto & Kr = make_field( + auto & Kr = muGrid::make_field( "Tangent Moduli reference", globalfields); // to be computed with tangent globalfields.initialise(cube, loc); static_assert( std::is_same::value, "oh oh"); static_assert( std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert( std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); using traits = MaterialMuSpectre_traits; { // block to contain not-constant gradient map typename traits::StressMap_t grad_map( globalfields["Transformation Gradient"]); for (auto F_ : grad_map) { F_.setRandom(); } grad_map[0] = grad_map[0].Identity(); // identifiable gradients for debug grad_map[1] = 1.2 * grad_map[1].Identity(); // ditto } // compute stresses using material mat.compute_stresses(globalfields["Transformation Gradient"], globalfields["Nominal Stress1"], Formulation::finite_strain); // compute stresses and tangent moduli using material BOOST_CHECK_THROW( mat.compute_stresses_tangent(globalfields["Transformation Gradient"], globalfields["Nominal Stress2"], globalfields["Nominal Stress2"], Formulation::finite_strain), std::runtime_error); mat.compute_stresses_tangent(globalfields["Transformation Gradient"], globalfields["Nominal Stress2"], globalfields["Tangent Moduli"], Formulation::finite_strain); typename traits::StrainMap_t Fmap(globalfields["Transformation Gradient"]); typename traits::StressMap_t Pmap_ref( globalfields["Nominal Stress reference"]); typename traits::TangentMap_t Kmap_ref( globalfields["Tangent Moduli reference"]); for (auto tup : akantu::zip(Fmap, Pmap_ref, Kmap_ref)) { auto F_ = std::get<0>(tup); auto P_ = std::get<1>(tup); auto K_ = std::get<2>(tup); - std::tie(P_, K_) = testGoodies::objective_hooke_explicit( - Fix::lambda, Fix::mu, F_); + std::tie(P_, K_) = + muGrid::testGoodies::objective_hooke_explicit(Fix::lambda, + Fix::mu, F_); } typename traits::StressMap_t Pmap_1(globalfields["Nominal Stress1"]); for (auto tup : akantu::zip(Pmap_ref, Pmap_1)) { auto P_r = std::get<0>(tup); auto P_1 = std::get<1>(tup); Real error = (P_r - P_1).norm(); BOOST_CHECK_LT(error, tol); } typename traits::StressMap_t Pmap_2(globalfields["Nominal Stress2"]); typename traits::TangentMap_t Kmap(globalfields["Tangent Moduli"]); for (auto tup : akantu::zip(Pmap_ref, Pmap_2, Kmap_ref, Kmap)) { auto P_r = std::get<0>(tup); auto P = std::get<1>(tup); Real error = (P_r - P).norm(); if (error >= tol) { std::cout << "error = " << error << " >= " << tol << " = tol" << std::endl; std::cout << "P(0) =" << std::endl << P_r << std::endl; std::cout << "P(1) =" << std::endl << P << std::endl; } BOOST_CHECK_LT(error, tol); auto K_r = std::get<2>(tup); auto K = std::get<3>(tup); error = (K_r - K).norm(); BOOST_CHECK_LT(error, tol); } } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_material_linear_elastic3.cc b/tests/test_material_linear_elastic3.cc index a25b0c8..6278737 100644 --- a/tests/test_material_linear_elastic3.cc +++ b/tests/test_material_linear_elastic3.cc @@ -1,92 +1,92 @@ /** * @file test_material_linear_elastic3.cc * * @author Richard Leute * * @date 21 Feb 2018 * * @brief description * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include "tests.hh" #include "materials/material_linear_elastic3.hh" #include "materials/materials_toolbox.hh" -#include "common/T4_map_proxy.hh" -#include "cmath" +#include +#include namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_linear_elastic_3); template struct MaterialFixture { using Material_t = Mat_t; Material_t mat; MaterialFixture() : mat("name") { mat.add_pixel({0}, Young, Poisson); } Real Young{10}; Real Poisson{0.3}; }; using mat_list = boost::mpl::list>, MaterialFixture>, MaterialFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constructor, Fix, mat_list, Fix){}; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_response, Fix, mat_list, Fix) { constexpr Dim_t Dim{Fix::Material_t::Parent::Parent::mdim()}; Eigen::Matrix E; E.setZero(); E(0, 0) = 0.001; E(1, 0) = E(0, 1) = 0.005; - using Hooke = - MatTB::Hooke, T4Mat>; + using Hooke = MatTB::Hooke, + muGrid::T4Mat>; Real lambda = Hooke::compute_lambda(Fix::Young, Fix::Poisson); Real mu = Hooke::compute_mu(Fix::Young, Fix::Poisson); auto C = Hooke::compute_C(lambda, mu); - T4MatMap Cmap{C.data()}; + muGrid::T4MatMap Cmap{C.data()}; Eigen::Matrix stress = Fix::mat.evaluate_stress(E, Cmap); Real sigma00 = lambda * E(0, 0) + 2 * mu * E(0, 0); Real sigma01 = 2 * mu * E(0, 1); Real sigma11 = lambda * E(0, 0); BOOST_CHECK_LT(std::abs(stress(0, 0) - sigma00), tol); BOOST_CHECK_LT(std::abs(stress(0, 1) - sigma01), tol); BOOST_CHECK_LT(std::abs(stress(1, 0) - sigma01), tol); BOOST_CHECK_LT(std::abs(stress(1, 1) - sigma11), tol); if (Dim == threeD) { for (int i = 0; i < Dim - 1; ++i) { BOOST_CHECK_LT(std::abs(stress(2, i)), tol); BOOST_CHECK_LT(std::abs(stress(i, 2)), tol); } BOOST_CHECK_LT(std::abs(stress(2, 2) - sigma11), tol); } }; BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_material_linear_elastic4.cc b/tests/test_material_linear_elastic4.cc index 5677597..d048898 100644 --- a/tests/test_material_linear_elastic4.cc +++ b/tests/test_material_linear_elastic4.cc @@ -1,97 +1,97 @@ /** * @file test_material_linear_elastic4.cc * * @author Richard Leute * * @date 27 Mar 2018 * * @brief description * * Copyright © 2018 Richard Leute * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include "tests.hh" #include "materials/material_linear_elastic4.hh" #include "materials/materials_toolbox.hh" -#include "common/T4_map_proxy.hh" -#include "cmath" +#include +#include namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_linear_elastic_4); template struct MaterialFixture { using Material_t = Mat_t; Material_t mat; MaterialFixture() : mat("name") { mat.add_pixel({0}, Youngs_modulus, Poisson_ratio); } Real Youngs_modulus{10}; Real Poisson_ratio{0.3}; }; using mat_list = boost::mpl::list>, MaterialFixture>, MaterialFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constructor, Fix, mat_list, Fix){}; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_response, Fix, mat_list, Fix) { constexpr Dim_t Dim{Fix::Material_t::Parent::Parent::mdim()}; Eigen::Matrix E; E.setZero(); E(0, 0) = 0.001; E(1, 0) = E(0, 1) = 0.005; - using Hooke = - MatTB::Hooke, T4Mat>; + using Hooke = MatTB::Hooke, + muGrid::T4Mat>; Real lambda = Hooke::compute_lambda(Fix::Youngs_modulus, Fix::Poisson_ratio); Real mu = Hooke::compute_mu(Fix::Youngs_modulus, Fix::Poisson_ratio); Eigen::Matrix stress = Fix::mat.evaluate_stress(E, lambda, mu); Real sigma00 = lambda * E(0, 0) + 2 * mu * E(0, 0); Real sigma01 = 2 * mu * E(0, 1); Real sigma11 = lambda * E(0, 0); BOOST_CHECK_LT(std::abs(stress(0, 0) - sigma00), tol); BOOST_CHECK_LT(std::abs(stress(0, 1) - sigma01), tol); BOOST_CHECK_LT(std::abs(stress(1, 0) - sigma01), tol); BOOST_CHECK_LT(std::abs(stress(1, 1) - sigma11), tol); if (Dim == threeD) { for (int i = 0; i < Dim - 1; ++i) { BOOST_CHECK_LT(std::abs(stress(2, i)), tol); BOOST_CHECK_LT(std::abs(stress(i, 2)), tol); } BOOST_CHECK_LT(std::abs(stress(2, 2) - sigma11), tol); } }; BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_material_linear_elastic_generic.cc b/tests/test_material_linear_elastic_generic.cc index 315ab1f..2a1b3dd 100644 --- a/tests/test_material_linear_elastic_generic.cc +++ b/tests/test_material_linear_elastic_generic.cc @@ -1,93 +1,93 @@ /** * @file test_material_linear_elastic_generic.cc * * @author Till Junge * * @date 21 Sep 2018 * * @brief test for the generic linear elastic law * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "tests.hh" #include "materials/material_linear_elastic_generic1.hh" #include "materials/materials_toolbox.hh" #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_linear_elastic_generic); template struct MatFixture { using Mat_t = MaterialLinearElasticGeneric1; using T2_t = Eigen::Matrix; - using T4_t = T4Mat; + using T4_t = muGrid::T4Mat; using V_t = Eigen::Matrix; constexpr static Real lambda{2}, mu{1.5}; constexpr static Real get_lambda() { return lambda; } constexpr static Real get_mu() { return mu; } constexpr static Real young{mu * (3 * lambda + 2 * mu) / (lambda + mu)}; constexpr static Real poisson{lambda / (2 * (lambda + mu))}; using Hooke = MatTB::Hooke; MatFixture() : C_voigt{get_C_voigt()}, mat("material", this->C_voigt) {} static V_t get_C_voigt() { V_t C{}; C.setZero(); C.template topLeftCorner().setConstant(get_lambda()); C.template topLeftCorner() += 2 * get_mu() * T2_t::Identity(); constexpr Dim_t Rest{vsize(Dim) - Dim}; using Rest_t = Eigen::Matrix; C.template bottomRightCorner() += get_mu() * Rest_t::Identity(); return C; } V_t C_voigt; Mat_t mat; }; using mats = boost::mpl::list, MatFixture>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(C_test, Fix, mats, Fix) { const auto ref_C{ Fix::Hooke::compute_C_T4(Fix::get_lambda(), Fix::get_mu())}; Real error{(ref_C - Fix::mat.get_C()).norm()}; BOOST_CHECK_LT(error, tol); if (not(error < tol)) { std::cout << "ref:" << std::endl << ref_C << std::endl; std::cout << "new:" << std::endl << Fix::mat.get_C() << std::endl; } } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_materials_toolbox.cc b/tests/test_materials_toolbox.cc index 729f0ff..2571d7d 100644 --- a/tests/test_materials_toolbox.cc +++ b/tests/test_materials_toolbox.cc @@ -1,371 +1,371 @@ /** * @file test_materials_toolbox.cc * * @author Till Junge * * @date 05 Nov 2017 * * @brief Tests for the materials toolbox * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include #include "tests.hh" #include "materials/materials_toolbox.hh" #include "materials/stress_transformations_default_case.hh" #include "materials/stress_transformations_PK1_impl.hh" #include "materials/stress_transformations_PK2_impl.hh" #include "materials/stress_transformations.hh" #include "tests/test_goodies.hh" #include #include -using namespace muGrid; - namespace muSpectre { BOOST_AUTO_TEST_SUITE(materials_toolbox) BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_strain_conversion, Fix, - testGoodies::dimlist, Fix) { + muGrid::testGoodies::dimlist, Fix) { constexpr Dim_t dim{Fix::dim}; using T2 = Eigen::Matrix; T2 F{(T2::Random() - .5 * T2::Ones()) / 20 + T2::Identity()}; // checking Green-Lagrange T2 Eref = .5 * (F.transpose() * F - T2::Identity()); T2 E_tb = MatTB::convert_strain( Eigen::Map>(F.data())); Real error = (Eref - E_tb).norm(); BOOST_CHECK_LT(error, tol); // checking Left Cauchy-Green Eref = F * F.transpose(); E_tb = MatTB::convert_strain(F); error = (Eref - E_tb).norm(); BOOST_CHECK_LT(error, tol); // checking Right Cauchy-Green Eref = F.transpose() * F; E_tb = MatTB::convert_strain(F); error = (Eref - E_tb).norm(); BOOST_CHECK_LT(error, tol); // checking Hencky (logarithmic) strain Eref = F.transpose() * F; Eigen::SelfAdjointEigenSolver EigSolv(Eref); Eref.setZero(); for (size_t i{0}; i < dim; ++i) { auto && vec = EigSolv.eigenvectors().col(i); auto && val = EigSolv.eigenvalues()(i); Eref += .5 * std::log(val) * vec * vec.transpose(); } E_tb = MatTB::convert_strain(F); error = (Eref - E_tb).norm(); BOOST_CHECK_LT(error, tol); auto F_tb = MatTB::convert_strain( F); error = (F - F_tb).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(dumb_tensor_mult_test, Fix, - testGoodies::dimlist, Fix) { + muGrid::testGoodies::dimlist, Fix) { constexpr Dim_t dim{Fix::dim}; - using T4 = T4Mat; + using T4 = muGrid::T4Mat; T4 A, B, R1, R2; A.setRandom(); B.setRandom(); R1 = A * B; R2.setZero(); for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t a = 0; a < dim; ++a) { for (Dim_t b = 0; b < dim; ++b) { for (Dim_t k = 0; k < dim; ++k) { for (Dim_t l = 0; l < dim; ++l) { - get(R2, i, j, k, l) += get(A, i, j, a, b) * get(B, a, b, k, l); + muGrid::get(R2, i, j, k, l) += + muGrid::get(A, i, j, a, b) * muGrid::get(B, a, b, k, l); } } } } } } auto error{(R1 - R2).norm()}; BOOST_CHECK_LT(error, tol); } - BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_PK1_stress, Fix, testGoodies::dimlist, - Fix) { + BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_PK1_stress, Fix, + muGrid::testGoodies::dimlist, Fix) { using Matrices::Tens2_t; using Matrices::Tens4_t; using Matrices::tensmult; constexpr Dim_t dim{Fix::dim}; using T2 = Eigen::Matrix; - using T4 = T4Mat; - testGoodies::RandRange rng; + using T4 = muGrid::T4Mat; + muGrid::testGoodies::RandRange rng; T2 F = T2::Identity() * 2; // F.setRandom(); T2 E_tb = MatTB::convert_strain( Eigen::Map>(F.data())); Real lambda = 3; // rng.randval(1, 2); Real mu = 4; // rng.randval(1,2); T4 J = Matrices::Itrac(); T2 I = Matrices::I2(); T4 I4 = Matrices::Isymm(); T4 C = lambda * J + 2 * mu * I4; T2 S = Matrices::tensmult(C, E_tb); T2 Sref = lambda * E_tb.trace() * I + 2 * mu * E_tb; auto error{(Sref - S).norm()}; BOOST_CHECK_LT(error, tol); T4 K = Matrices::outer_under(I, S) + Matrices::outer_under(F, I) * C * Matrices::outer_under(F.transpose(), I); // See Curnier, 2000, "Méthodes numériques en mécanique des solides", p 252 T4 Kref; Real Fkrkr = (F.array() * F.array()).sum(); T2 Fkmkn = F.transpose() * F; T2 Fisjs = F * F.transpose(); Kref.setZero(); for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t m = 0; m < dim; ++m) { for (Dim_t n = 0; n < dim; ++n) { - get(Kref, i, m, j, n) = + muGrid::get(Kref, i, m, j, n) = (lambda * ((Fkrkr - dim) / 2 * I(i, j) * I(m, n) + F(i, m) * F(j, n)) + mu * (I(i, j) * Fkmkn(m, n) + Fisjs(i, j) * I(m, n) - I(i, j) * I(m, n) + F(i, n) * F(j, m))); } } } } error = (Kref - K).norm(); BOOST_CHECK_LT(error, tol); T2 P = MatTB::PK1_stress( F, S); T2 Pref = F * S; error = (P - Pref).norm(); BOOST_CHECK_LT(error, tol); auto && stress_tgt = MatTB::PK1_stress( F, S, C); T2 P_t = std::move(std::get<0>(stress_tgt)); T4 K_t = std::move(std::get<1>(stress_tgt)); error = (P_t - Pref).norm(); BOOST_CHECK_LT(error, tol); error = (K_t - Kref).norm(); BOOST_CHECK_LT(error, tol); auto && stress_tgt_trivial = MatTB::PK1_stress(F, P, K); T2 P_u = std::move(std::get<0>(stress_tgt_trivial)); T4 K_u = std::move(std::get<1>(stress_tgt_trivial)); error = (P_u - Pref).norm(); BOOST_CHECK_LT(error, tol); error = (K_u - Kref).norm(); BOOST_CHECK_LT(error, tol); T2 P_g; T4 K_g; - std::tie(P_g, K_g) = testGoodies::objective_hooke_explicit(lambda, mu, F); + std::tie(P_g, K_g) = + muGrid::testGoodies::objective_hooke_explicit(lambda, mu, F); error = (P_g - Pref).norm(); BOOST_CHECK_LT(error, tol); error = (K_g - Kref).norm(); BOOST_CHECK_LT(error, tol); } BOOST_AUTO_TEST_CASE(elastic_modulus_conversions) { // define original input constexpr Real E{123.456}; constexpr Real nu{.3}; using MatTB::convert_elastic_modulus; // derived values constexpr Real K{ convert_elastic_modulus(E, nu)}; constexpr Real lambda{ convert_elastic_modulus(E, nu)}; constexpr Real mu{ convert_elastic_modulus(E, nu)}; // recover original inputs Real comp = convert_elastic_modulus(K, mu); Real err = E - comp; BOOST_CHECK_LT(err, tol); comp = convert_elastic_modulus(K, mu); err = nu - comp; BOOST_CHECK_LT(err, tol); comp = convert_elastic_modulus(lambda, mu); err = E - comp; BOOST_CHECK_LT(err, tol); // check inversion resistance Real compA = convert_elastic_modulus(K, mu); Real compB = convert_elastic_modulus(mu, K); BOOST_CHECK_EQUAL(compA, compB); // check trivial self-returning comp = convert_elastic_modulus(K, mu); BOOST_CHECK_EQUAL(K, comp); comp = convert_elastic_modulus(K, mu); BOOST_CHECK_EQUAL(mu, comp); // check alternative calculation of computed values comp = convert_elastic_modulus( K, mu); // alternative for "Shear" BOOST_CHECK_LE(std::abs((comp - lambda) / lambda), tol); } template struct FiniteDifferencesHolder { constexpr static FiniteDiff value{FinDiff}; }; using FinDiffList = boost::mpl::list, FiniteDifferencesHolder, FiniteDifferencesHolder>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(numerical_tangent_test, Fix, FinDiffList, Fix) { constexpr Dim_t Dim{twoD}; - using T4_t = T4Mat; + using T4_t = muGrid::T4Mat; using T2_t = Eigen::Matrix; bool verbose{false}; T4_t Q{}; Q << 1., 2., 0., 0., 0., 1.66666667, 0., 0., 0., 0., 2.33333333, 0., 0., 0., 0., 3.; if (verbose) { std::cout << Q << std::endl << std::endl; } T2_t B{}; B << 2., 3.33333333, 2.66666667, 4.; if (verbose) { std::cout << B << std::endl << std::endl; } auto fun = [&](const T2_t & x) -> T2_t { using cmap_t = Eigen::Map>; using map_t = Eigen::Map>; cmap_t x_vec{x.data()}; T2_t ret_val{}; map_t(ret_val.data()) = Q * x_vec + cmap_t(B.data()); return ret_val; }; T2_t temp_res = fun(T2_t::Ones()); if (verbose) { std::cout << temp_res << std::endl << std::endl; } T4_t numerical_tangent{MatTB::compute_numerical_tangent( fun, T2_t::Ones(), 1e-2)}; if (verbose) { std::cout << numerical_tangent << std::endl << std::endl; } Real error{(numerical_tangent - Q).norm() / Q.norm()}; BOOST_CHECK_LT(error, tol); if (not(error < tol)) { switch (Fix::value) { case FiniteDiff::backward: { std::cout << "backward difference: " << std::endl; break; } case FiniteDiff::forward: { std::cout << "forward difference: " << std::endl; break; } case FiniteDiff::centred: { std::cout << "centered difference: " << std::endl; break; } } std::cout << "error = " << error << std::endl; std::cout << "numerical tangent:\n" << numerical_tangent << std::endl; std::cout << "reference:\n" << Q << std::endl; } } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_projection.hh b/tests/test_projection.hh index 3e8b9bd..0020ff4 100644 --- a/tests/test_projection.hh +++ b/tests/test_projection.hh @@ -1,105 +1,106 @@ /** * @file test_projection.hh * * @author Till Junge * * @date 16 Jan 2018 * * @brief common declarations for testing both the small and finite strain * projection operators * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "tests.hh" -#include "projection/fftw_engine.hh" +#include #include #include +#include #ifndef TESTS_TEST_PROJECTION_HH_ #define TESTS_TEST_PROJECTION_HH_ namespace muSpectre { /* ---------------------------------------------------------------------- */ template struct Sizes {}; template <> struct Sizes { constexpr static Ccoord_t get_resolution() { return Ccoord_t{3, 5}; } constexpr static Rcoord_t get_lengths() { return Rcoord_t{3.4, 5.8}; } }; template <> struct Sizes { constexpr static Ccoord_t get_resolution() { return Ccoord_t{3, 5, 7}; } constexpr static Rcoord_t get_lengths() { return Rcoord_t{3.4, 5.8, 6.7}; } }; template struct Squares {}; template <> struct Squares { constexpr static Ccoord_t get_resolution() { return Ccoord_t{5, 5}; } constexpr static Rcoord_t get_lengths() { return Rcoord_t{5, 5}; } }; template <> struct Squares { constexpr static Ccoord_t get_resolution() { return Ccoord_t{7, 7, 7}; } constexpr static Rcoord_t get_lengths() { return Rcoord_t{7, 7, 7}; } }; /* ---------------------------------------------------------------------- */ template struct ProjectionFixture { - using Engine = FFTWEngine; + using Engine = muFFT::FFTWEngine; using Parent = Proj; constexpr static Dim_t sdim{DimS}; constexpr static Dim_t mdim{DimM}; ProjectionFixture() : projector(std::make_unique(SizeGiver::get_resolution(), mdim * mdim), SizeGiver::get_lengths()) {} Parent projector; }; } // namespace muSpectre #endif // TESTS_TEST_PROJECTION_HH_ diff --git a/tests/test_projection_finite.cc b/tests/test_projection_finite.cc index 6064e66..9109ea4 100644 --- a/tests/test_projection_finite.cc +++ b/tests/test_projection_finite.cc @@ -1,150 +1,151 @@ /** * @file test_projection_finite.cc * * @author Till Junge * * @date 07 Dec 2017 * * @brief tests for standard finite strain projection operator * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "projection/projection_finite_strain.hh" #include "projection/projection_finite_strain_fast.hh" -#include "projection/fft_utils.hh" +#include #include "test_projection.hh" #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(projection_finite_strain); /* ---------------------------------------------------------------------- */ using fixlist = boost::mpl::list< ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrainFast>, ProjectionFixture, ProjectionFiniteStrainFast>, ProjectionFixture, ProjectionFiniteStrainFast>, ProjectionFixture, ProjectionFiniteStrainFast>>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { - BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); + BOOST_CHECK_NO_THROW( + fix::projector.initialise(muFFT::FFT_PlanFlags::estimate)); } /* ---------------------------------------------------------------------- */ BOOST_AUTO_TEST_CASE(even_grid_test) { - using Engine = FFTWEngine; + using Engine = muFFT::FFTWEngine; using proj = ProjectionFiniteStrainFast; auto engine = std::make_unique(Ccoord_t{2, 2}, 2 * 2); BOOST_CHECK_THROW(proj(std::move(engine), Rcoord_t{4.3, 4.3}), std::runtime_error); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test, fix, fixlist, fix) { // create a gradient field with a zero mean gradient and verify // that the projection preserves it constexpr Dim_t dim{fix::sdim}, sdim{fix::sdim}, mdim{fix::mdim}; static_assert( dim == fix::mdim, "These tests assume that the material and spatial dimension are " "identical"); - using Fields = GlobalFieldCollection; - using FieldT = TensorField; - using FieldMap = MatrixFieldMap; + using Fields = muGrid::GlobalFieldCollection; + using FieldT = muGrid::TensorField; + using FieldMap = muGrid::MatrixFieldMap; using Vector = Eigen::Matrix; Fields fields{}; - FieldT & f_grad{make_field("gradient", fields)}; - FieldT & f_var{make_field("working field", fields)}; + FieldT & f_grad{muGrid::make_field("gradient", fields)}; + FieldT & f_var{muGrid::make_field("working field", fields)}; FieldMap grad(f_grad); FieldMap var(f_var); fields.initialise(fix::projector.get_subdomain_resolutions(), fix::projector.get_subdomain_locations()); - FFT_freqs freqs{fix::projector.get_domain_resolutions(), - fix::projector.get_domain_lengths()}; + muFFT::FFT_freqs freqs{fix::projector.get_domain_resolutions(), + fix::projector.get_domain_lengths()}; Vector k; for (Dim_t i = 0; i < dim; ++i) { // the wave vector has to be such that it leads to an integer // number of periods in each length of the domain - k(i) = (i + 1) * 2 * pi / fix::projector.get_domain_lengths()[i]; + k(i) = (i + 1) * 2 * muGrid::pi / fix::projector.get_domain_lengths()[i]; } + using muGrid::operator/; for (auto && tup : akantu::zip(fields, grad, var)) { auto & ccoord = std::get<0>(tup); auto & g = std::get<1>(tup); auto & v = std::get<2>(tup); - Vector vec = CcoordOps::get_vector( + Vector vec = muGrid::CcoordOps::get_vector( ccoord, fix::projector.get_domain_lengths() / fix::projector.get_domain_resolutions()); g.row(0) = k.transpose() * cos(k.dot(vec)); v.row(0) = g.row(0); } - fix::projector.initialise(FFT_PlanFlags::estimate); + fix::projector.initialise(muFFT::FFT_PlanFlags::estimate); fix::projector.apply_projection(f_var); for (auto && tup : akantu::zip(fields, grad, var)) { auto & ccoord = std::get<0>(tup); auto & g = std::get<1>(tup); auto & v = std::get<2>(tup); - Vector vec = CcoordOps::get_vector( + Vector vec = muGrid::CcoordOps::get_vector( ccoord, fix::projector.get_domain_lengths() / fix::projector.get_domain_resolutions()); Real error = (g - v).norm(); BOOST_CHECK_LT(error, tol); if (error >= tol) { std::cout << std::endl << "grad_ref :" << std::endl << g << std::endl; std::cout << std::endl << "grad_proj :" << std::endl << v << std::endl; - std::cout << std::endl - << "ccoord :" << std::endl - << ccoord << std::endl; + std::cout << std::endl << "ccoord :" << std::endl; + muGrid::operator<<(std::cout, ccoord) << std::endl; std::cout << std::endl << "vector :" << std::endl << vec.transpose() << std::endl; } } } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_projection_small.cc b/tests/test_projection_small.cc index acc9ddd..0e57519 100644 --- a/tests/test_projection_small.cc +++ b/tests/test_projection_small.cc @@ -1,143 +1,145 @@ /** * @file test_projection_small.cc * * @author Till Junge * * @date 16 Jan 2018 * * @brief tests for standard small strain projection operator * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ -#include "projection/projection_small_strain.hh" #include "test_projection.hh" -#include "projection/fft_utils.hh" +#include "projection/projection_small_strain.hh" +#include #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(projection_small_strain); using fixlist = boost::mpl::list< ProjectionFixture, ProjectionSmallStrain>, ProjectionFixture, ProjectionSmallStrain>, ProjectionFixture, ProjectionSmallStrain>, ProjectionFixture, ProjectionSmallStrain>>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { - BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); + BOOST_CHECK_NO_THROW( + fix::projector.initialise(muFFT::FFT_PlanFlags::estimate)); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test, fix, fixlist, fix) { // create a gradient field with a zero mean gradient and verify // that the projection preserves it constexpr Dim_t dim{fix::sdim}, sdim{fix::sdim}, mdim{fix::mdim}; static_assert( dim == fix::mdim, "These tests assume that the material and spatial dimension are " "identical"); - using Fields = GlobalFieldCollection; - using FieldT = TensorField; - using FieldMap = MatrixFieldMap; + using Fields = muGrid::GlobalFieldCollection; + using FieldT = muGrid::TensorField; + using FieldMap = muGrid::MatrixFieldMap; using Vector = Eigen::Matrix; Fields fields{}; - FieldT & f_grad{make_field("strain", fields)}; - FieldT & f_var{make_field("working field", fields)}; + FieldT & f_grad{muGrid::make_field("strain", fields)}; + FieldT & f_var{muGrid::make_field("working field", fields)}; FieldMap grad(f_grad); FieldMap var(f_var); fields.initialise(fix::projector.get_subdomain_resolutions(), fix::projector.get_subdomain_locations()); - FFT_freqs freqs{fix::projector.get_domain_resolutions(), - fix::projector.get_domain_lengths()}; + muFFT::FFT_freqs freqs{fix::projector.get_domain_resolutions(), + fix::projector.get_domain_lengths()}; Vector k; for (Dim_t i = 0; i < dim; ++i) { // the wave vector has to be such that it leads to an integer // number of periods in each length of the domain - k(i) = (i + 1) * 2 * pi / fix::projector.get_domain_lengths()[i]; + k(i) = (i + 1) * 2 * muGrid::pi / fix::projector.get_domain_lengths()[i]; } + using muGrid::operator/; for (auto && tup : akantu::zip(fields, grad, var)) { auto & ccoord = std::get<0>(tup); auto & g = std::get<1>(tup); auto & v = std::get<2>(tup); - Vector vec = CcoordOps::get_vector( + Vector vec = muGrid::CcoordOps::get_vector( ccoord, fix::projector.get_domain_lengths() / fix::projector.get_domain_resolutions()); g.row(0) << k.transpose() * cos(k.dot(vec)); // We need to add I to the term, because this field has a net // zero gradient, which leads to a net -I strain g = 0.5 * ((g - g.Identity()).transpose() + (g - g.Identity())).eval() + g.Identity(); v = g; } - fix::projector.initialise(FFT_PlanFlags::estimate); + fix::projector.initialise(muFFT::FFT_PlanFlags::estimate); fix::projector.apply_projection(f_var); + using muGrid::operator/; constexpr bool verbose{false}; for (auto && tup : akantu::zip(fields, grad, var)) { auto & ccoord = std::get<0>(tup); auto & g = std::get<1>(tup); auto & v = std::get<2>(tup); - Vector vec = CcoordOps::get_vector( + Vector vec = muGrid::CcoordOps::get_vector( ccoord, fix::projector.get_domain_lengths() / fix::projector.get_domain_resolutions()); Real error = (g - v).norm(); BOOST_CHECK_LT(error, tol); if ((error >= tol) || verbose) { std::cout << std::endl << "grad_ref :" << std::endl << g << std::endl; std::cout << std::endl << "grad_proj :" << std::endl << v << std::endl; - std::cout << std::endl - << "ccoord :" << std::endl - << ccoord << std::endl; + std::cout << std::endl << "ccoord :" << std::endl; + muGrid::operator<<(std::cout, ccoord) << std::endl; std::cout << std::endl << "vector :" << std::endl << vec.transpose() << std::endl; std::cout << "means:" << std::endl << ":" << std::endl << grad.mean() << std::endl << ":" << std::endl << var.mean(); } } } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre diff --git a/tests/test_solver_newton_cg.cc b/tests/test_solver_newton_cg.cc index f0228a5..fca2a6a 100644 --- a/tests/test_solver_newton_cg.cc +++ b/tests/test_solver_newton_cg.cc @@ -1,481 +1,483 @@ /** * @file test_solver_newton_cg.cc * * @author Till Junge * * @date 20 Dec 2017 * * @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 Lesser 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #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 "projection/fftw_engine.hh" #include "projection/projection_finite_strain_fast.hh" #include "materials/material_linear_elastic1.hh" -#include "common/iterators.hh" -#include "common/ccoord_operations.hh" #include "cell/cell_factory.hh" +#include +#include +#include + #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(newton_cg_tests); BOOST_AUTO_TEST_CASE(manual_construction_test) { // constexpr Dim_t dim{twoD}; constexpr Dim_t dim{threeD}; // constexpr Ccoord_t resolutions{3, 3}; // constexpr Rcoord_t lengths{2.3, 2.7}; constexpr Ccoord_t resolutions{5, 5, 5}; constexpr Rcoord_t lengths{5, 5, 5}; - auto fft_ptr{std::make_unique>(resolutions, ipow(dim, 2))}; + auto fft_ptr{std::make_unique>( + resolutions, muGrid::ipow(dim, 2))}; auto proj_ptr{std::make_unique>( std::move(fft_ptr), lengths)}; CellBase sys(std::move(proj_ptr)); using Mat_t = MaterialLinearElastic1; // const Real Young{210e9}, Poisson{.33}; const Real Young{1.0030648180242636}, Poisson{0.29930675909878679}; // const Real lambda{Young*Poisson/((1+Poisson)*(1-2*Poisson))}; // const Real mu{Young/(2*(1+Poisson))}; auto & Material_hard = Mat_t::make(sys, "hard", 10 * Young, Poisson); auto & Material_soft = Mat_t::make(sys, "soft", Young, Poisson); for (auto && tup : akantu::enumerate(sys)) { auto && pixel = std::get<1>(tup); if (std::get<0>(tup) == 0) { Material_hard.add_pixel(pixel); } else { Material_soft.add_pixel(pixel); } } sys.initialise(); Grad_t delF0; delF0 << 0, 1., 0, 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 Uint maxiter{muGrid::CcoordOps::get_size(resolutions) * + muGrid::ipow(dim, secondOrder) * 10}; constexpr bool verbose{false}; GradIncrements grads; grads.push_back(delF0); DeprecatedSolverCG cg{sys, cg_tol, maxiter, static_cast(verbose)}; Eigen::ArrayXXd res1{ deprecated_de_geus(sys, grads, cg, newton_tol, verbose)[0].grad}; DeprecatedSolverCG cg2{sys, cg_tol, maxiter, static_cast(verbose)}; Eigen::ArrayXXd res2{ deprecated_newton_cg(sys, grads, cg2, newton_tol, verbose)[0].grad}; BOOST_CHECK_LE(abs(res1 - res2).mean(), cg_tol); } BOOST_AUTO_TEST_CASE(small_strain_patch_test) { constexpr Dim_t dim{twoD}; using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; - constexpr Ccoord resolutions{CcoordOps::get_cube(3)}; - constexpr Rcoord lengths{CcoordOps::get_cube(1.)}; + constexpr Ccoord resolutions{muGrid::CcoordOps::get_cube(3)}; + constexpr Rcoord lengths{muGrid::CcoordOps::get_cube(1.)}; constexpr Formulation form{Formulation::small_strain}; // number of layers in the hard material constexpr Uint nb_lays{1}; constexpr Real contrast{2}; static_assert(nb_lays < resolutions[0], "the number or layers in the hard material must be smaller " "than the total number of layers in dimension 0"); auto sys{make_cell(resolutions, lengths, form)}; using Mat_t = MaterialLinearElastic1; constexpr Real Young{2.}, Poisson{.33}; auto material_hard{ std::make_unique("hard", contrast * Young, Poisson)}; auto material_soft{std::make_unique("soft", Young, Poisson)}; for (const auto & pixel : sys) { if (pixel[0] < Dim_t(nb_lays)) { material_hard->add_pixel(pixel); } else { material_soft->add_pixel(pixel); } } sys.add_material(std::move(material_hard)); sys.add_material(std::move(material_soft)); sys.initialise(); Grad_t delEps0{Grad_t::Zero()}; constexpr Real eps0 = 1.; // delEps0(0, 1) = delEps0(1, 0) = eps0; delEps0(0, 0) = eps0; constexpr Real cg_tol{1e-8}, newton_tol{1e-5}, equil_tol{1e-10}; constexpr Uint maxiter{dim * 10}; constexpr Dim_t verbose{0}; DeprecatedSolverCGEigen cg{sys, cg_tol, maxiter, static_cast(verbose)}; auto result = deprecated_newton_cg( sys, delEps0, cg, newton_tol, // de_geus(sys, delEps0, cg, newton_tol, equil_tol, verbose); if (verbose) { std::cout << "result:" << std::endl << result.grad << std::endl; std::cout << "mean strain = " << std::endl << sys.get_strain().get_map().mean() << std::endl; } /** * verification of resultant strains: subscript ₕ for hard and ₛ * for soft, Nₕ is nb_lays and Nₜₒₜ is resolutions, k is contrast * * Δl = εl = Δlₕ + Δlₛ = εₕlₕ+εₛlₛ * => ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ * * σ is constant across all layers * σₕ = σₛ * => Eₕ εₕ = Eₛ εₛ * => εₕ = 1/k εₛ * => ε / (1/k Nₕ/Nₜₒₜ + (Nₜₒₜ-Nₕ)/Nₜₒₜ) = εₛ */ constexpr Real factor{1 / contrast * Real(nb_lays) / resolutions[0] + 1. - nb_lays / Real(resolutions[0])}; constexpr Real eps_soft{eps0 / factor}; constexpr Real eps_hard{eps_soft / contrast}; if (verbose) { std::cout << "εₕ = " << eps_hard << ", εₛ = " << eps_soft << std::endl; std::cout << "ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ" << std::endl; } Grad_t Eps_hard; Eps_hard << eps_hard, 0, 0, 0; Grad_t Eps_soft; Eps_soft << eps_soft, 0, 0, 0; // verify uniaxial tension patch test for (const auto & pixel : sys) { if (pixel[0] < Dim_t(nb_lays)) { BOOST_CHECK_LE((Eps_hard - sys.get_strain().get_map()[pixel]).norm(), tol); } else { BOOST_CHECK_LE((Eps_soft - sys.get_strain().get_map()[pixel]).norm(), tol); } } delEps0 = Grad_t::Zero(); delEps0(0, 1) = delEps0(1, 0) = eps0; DeprecatedSolverCG cg2{sys, cg_tol, maxiter, static_cast(verbose)}; result = deprecated_newton_cg(sys, delEps0, cg2, newton_tol, equil_tol, verbose); Eps_hard << 0, eps_hard, eps_hard, 0; Eps_soft << 0, eps_soft, eps_soft, 0; // verify pure shear patch test for (const auto & pixel : sys) { if (pixel[0] < Dim_t(nb_lays)) { BOOST_CHECK_LE((Eps_hard - sys.get_strain().get_map()[pixel]).norm(), tol); } else { BOOST_CHECK_LE((Eps_soft - sys.get_strain().get_map()[pixel]).norm(), tol); } } } template struct SolverFixture { using type = SolverType; }; using SolverList = boost::mpl::list< SolverFixture, SolverFixture, SolverFixture, SolverFixture, SolverFixture, SolverFixture>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(small_strain_patch_dynamic_solver, Fix, SolverList, Fix) { // BOOST_AUTO_TEST_CASE(small_strain_patch_test_dynamic_solver) { constexpr Dim_t dim{twoD}; using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; - constexpr Ccoord resolutions{CcoordOps::get_cube(3)}; - constexpr Rcoord lengths{CcoordOps::get_cube(1.)}; + constexpr Ccoord resolutions{muGrid::CcoordOps::get_cube(3)}; + constexpr Rcoord lengths{muGrid::CcoordOps::get_cube(1.)}; constexpr Formulation form{Formulation::small_strain}; // number of layers in the hard material constexpr Uint nb_lays{1}; constexpr Real contrast{2}; static_assert(nb_lays < resolutions[0], "the number or layers in the hard material must be smaller " "than the total number of layers in dimension 0"); auto sys{make_cell(resolutions, lengths, form)}; using Mat_t = MaterialLinearElastic1; constexpr Real Young{2.}, Poisson{.33}; auto material_hard{ std::make_unique("hard", contrast * Young, Poisson)}; auto material_soft{std::make_unique("soft", Young, Poisson)}; for (const auto & pixel : sys) { if (pixel[0] < Dim_t(nb_lays)) { material_hard->add_pixel(pixel); } else { material_soft->add_pixel(pixel); } } sys.add_material(std::move(material_hard)); sys.add_material(std::move(material_soft)); sys.initialise(); Grad_t delEps0{Grad_t::Zero()}; constexpr Real eps0 = 1.; // delEps0(0, 1) = delEps0(1, 0) = eps0; delEps0(0, 0) = eps0; constexpr Real cg_tol{1e-8}, newton_tol{1e-5}, equil_tol{1e-10}; constexpr Uint maxiter{dim * 10}; constexpr Dim_t verbose{0}; using Solver_t = typename Fix::type; Solver_t cg{sys, cg_tol, maxiter, static_cast(verbose)}; auto result = newton_cg(sys, delEps0, cg, newton_tol, equil_tol, verbose); if (verbose) { std::cout << "result:" << std::endl << result.grad << std::endl; std::cout << "mean strain = " << std::endl << sys.get_strain().get_map().mean() << std::endl; } /** * verification of resultant strains: subscript ₕ for hard and ₛ * for soft, Nₕ is nb_lays and Nₜₒₜ is resolutions, k is contrast * * Δl = εl = Δlₕ + Δlₛ = εₕlₕ+εₛlₛ * => ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ * * σ is constant across all layers * σₕ = σₛ * => Eₕ εₕ = Eₛ εₛ * => εₕ = 1/k εₛ * => ε / (1/k Nₕ/Nₜₒₜ + (Nₜₒₜ-Nₕ)/Nₜₒₜ) = εₛ */ constexpr Real factor{1 / contrast * Real(nb_lays) / resolutions[0] + 1. - nb_lays / Real(resolutions[0])}; constexpr Real eps_soft{eps0 / factor}; constexpr Real eps_hard{eps_soft / contrast}; if (verbose) { std::cout << "εₕ = " << eps_hard << ", εₛ = " << eps_soft << std::endl; std::cout << "ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ" << std::endl; } Grad_t Eps_hard; Eps_hard << eps_hard, 0, 0, 0; Grad_t Eps_soft; Eps_soft << eps_soft, 0, 0, 0; // verify uniaxial tension patch test for (const auto & pixel : sys) { if (pixel[0] < Dim_t(nb_lays)) { BOOST_CHECK_LE((Eps_hard - sys.get_strain().get_map()[pixel]).norm(), tol); } else { BOOST_CHECK_LE((Eps_soft - sys.get_strain().get_map()[pixel]).norm(), tol); } } delEps0 = Grad_t::Zero(); delEps0(0, 1) = delEps0(1, 0) = eps0; Solver_t cg2{sys, cg_tol, maxiter, static_cast(verbose)}; result = de_geus(sys, delEps0, cg2, newton_tol, equil_tol, verbose); Eps_hard << 0, eps_hard, eps_hard, 0; Eps_soft << 0, eps_soft, eps_soft, 0; // verify pure shear patch test for (const auto & pixel : sys) { if (pixel[0] < Dim_t(nb_lays)) { BOOST_CHECK_LE((Eps_hard - sys.get_strain().get_map()[pixel]).norm(), tol); } else { BOOST_CHECK_LE((Eps_soft - sys.get_strain().get_map()[pixel]).norm(), tol); } } } BOOST_AUTO_TEST_CASE(small_strain_patch_test_new_interface_manual) { constexpr Dim_t dim{twoD}; using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; - constexpr Ccoord resolutions{CcoordOps::get_cube(3)}; - constexpr Rcoord lengths{CcoordOps::get_cube(1.)}; + constexpr Ccoord resolutions{muGrid::CcoordOps::get_cube(3)}; + constexpr Rcoord lengths{muGrid::CcoordOps::get_cube(1.)}; constexpr Formulation form{Formulation::small_strain}; // number of layers in the hard material constexpr Uint nb_lays{1}; constexpr Real contrast{2}; static_assert(nb_lays < resolutions[0], "the number or layers in the hard material must be smaller " "than the total number of layers in dimension 0"); auto sys{make_cell(resolutions, lengths, form)}; using Mat_t = MaterialLinearElastic1; constexpr Real Young{2.}, Poisson{.33}; auto material_hard{ std::make_unique("hard", contrast * Young, Poisson)}; auto material_soft{std::make_unique("soft", Young, Poisson)}; for (const auto & pixel : sys) { if (pixel[0] < Dim_t(nb_lays)) { material_hard->add_pixel(pixel); } else { material_soft->add_pixel(pixel); } } sys.add_material(std::move(material_hard)); sys.add_material(std::move(material_soft)); Grad_t delEps0{Grad_t::Zero()}; constexpr Real eps0 = 1.; // delEps0(0, 1) = delEps0(1, 0) = eps0; delEps0(0, 0) = eps0; constexpr Real cg_tol{1e-8}; constexpr Uint maxiter{dim * 10}; constexpr Dim_t verbose{0}; DeprecatedSolverCGEigen cg{sys, cg_tol, maxiter, static_cast(verbose)}; auto F = sys.get_strain_vector(); F.setZero(); sys.evaluate_stress_tangent(); Eigen::VectorXd DelF(sys.get_nb_dof()); - using RMap_t = RawFieldMap>>; + using RMap_t = muGrid::RawFieldMap>>; for (auto tmp : RMap_t(DelF)) { tmp = delEps0; } Eigen::VectorXd rhs = -sys.evaluate_projected_directional_stiffness(DelF); F += DelF; DelF.setZero(); cg.initialise(); Eigen::Map(DelF.data(), DelF.size()) = cg.solve(rhs, DelF); F += DelF; if (verbose) { std::cout << "result:" << std::endl << F << std::endl; std::cout << "mean strain = " << std::endl << sys.get_strain().get_map().mean() << std::endl; } /** * verification of resultant strains: subscript ₕ for hard and ₛ * for soft, Nₕ is nb_lays and Nₜₒₜ is resolutions, k is contrast * * Δl = εl = Δlₕ + Δlₛ = εₕlₕ+εₛlₛ * => ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ * * σ is constant across all layers * σₕ = σₛ * => Eₕ εₕ = Eₛ εₛ * => εₕ = 1/k εₛ * => ε / (1/k Nₕ/Nₜₒₜ + (Nₜₒₜ-Nₕ)/Nₜₒₜ) = εₛ */ constexpr Real factor{1 / contrast * Real(nb_lays) / resolutions[0] + 1. - nb_lays / Real(resolutions[0])}; constexpr Real eps_soft{eps0 / factor}; constexpr Real eps_hard{eps_soft / contrast}; if (verbose) { std::cout << "εₕ = " << eps_hard << ", εₛ = " << eps_soft << std::endl; std::cout << "ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ" << std::endl; } Grad_t Eps_hard; Eps_hard << eps_hard, 0, 0, 0; Grad_t Eps_soft; Eps_soft << eps_soft, 0, 0, 0; // verify uniaxial tension patch test for (const auto & pixel : sys) { if (pixel[0] < Dim_t(nb_lays)) { BOOST_CHECK_LE((Eps_hard - sys.get_strain().get_map()[pixel]).norm(), tol); } else { BOOST_CHECK_LE((Eps_soft - sys.get_strain().get_map()[pixel]).norm(), tol); } } delEps0.setZero(); delEps0(0, 1) = delEps0(1, 0) = eps0; DeprecatedSolverCG cg2{sys, cg_tol, maxiter, static_cast(verbose)}; F.setZero(); sys.evaluate_stress_tangent(); for (auto tmp : RMap_t(DelF)) { tmp = delEps0; } rhs = -sys.evaluate_projected_directional_stiffness(DelF); F += DelF; DelF.setZero(); cg2.initialise(); DelF = cg2.solve(rhs, DelF); F += DelF; Eps_hard << 0, eps_hard, eps_hard, 0; Eps_soft << 0, eps_soft, eps_soft, 0; // verify pure shear patch test for (const auto & pixel : sys) { if (pixel[0] < Dim_t(nb_lays)) { BOOST_CHECK_LE((Eps_hard - sys.get_strain().get_map()[pixel]).norm(), tol); } else { BOOST_CHECK_LE((Eps_soft - sys.get_strain().get_map()[pixel]).norm(), tol); } } } BOOST_AUTO_TEST_SUITE_END(); } // namespace muSpectre