diff --git a/bin/demonstrator1.cc b/bin/demonstrator1.cc index 5fd2edc..62df04e 100644 --- a/bin/demonstrator1.cc +++ b/bin/demonstrator1.cc @@ -1,151 +1,142 @@ /** * file demonstrator1.cc * * @author Till Junge * * @date 03 Jan 2018 * * @brief larger problem to show off * * @section LICENCE * * 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 #include #include #include "external/cxxopts.hpp" #include "common/common.hh" #include "common/ccoord_operations.hh" #include "system/system_factory.hh" #include "materials/material_hyper_elastic1.hh" #include "solver/solvers.hh" #include "solver/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 muSpectre; int main(int argc, char *argv[]) { - - std::cout - << "µSpectre demonstrator1 Copyright © 2018 Till Junge " - << std::endl - << "This program comes with ABSOLUTELY NO WARRANTY." - << std::endl - << "This is free software, and you are welcome to redistribute it" - << std::endl - << "under certain conditions, see the license file." - << std::endl << std::endl; + 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 system{make_system(resolutions, lengths)}; constexpr Real E{1.0030648180242636}; constexpr Real nu{0.29930675909878679}; using Material_t = MaterialHyperElastic1; auto Material_soft{std::make_unique("soft", E, nu)}; auto Material_hard{std::make_unique("hard", 10*E, nu)}; int counter{0}; for (const auto && pixel:system) { 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 " << system.size() << " are in the hard material" << std::endl; system.add_material(std::move(Material_soft)); system.add_material(std::move(Material_hard)); system.initialise(FFT_PlanFlags::measure); constexpr Real newton_tol{1e-4}; constexpr Real cg_tol{1e-7}; const size_t 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}; de_geus(system, grads, form, cg_tol, newton_tol, maxiter, 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/src/common/common.cc b/src/common/common.cc index 89f3fcb..4673240 100644 --- a/src/common/common.cc +++ b/src/common/common.cc @@ -1,85 +1,100 @@ /** * file common.cc * * @author Till Junge * * @date 15 Nov 2017 * * @brief Implementation for common functions * * @section LICENCE * * 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 GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include namespace muSpectre { /* ---------------------------------------------------------------------- */ std::ostream & operator<<(std::ostream & os, Formulation f) { switch (f) { case Formulation::small_strain: {os << "small_strain"; break;} case Formulation::finite_strain: {os << "finite_strain"; break;} default: throw std::runtime_error ("unknown formulation."); break; } return os; } /* ---------------------------------------------------------------------- */ std::ostream & operator<<(std::ostream & os, StressMeasure s) { switch (s) { case StressMeasure::Cauchy: {os << "Cauchy"; break;} case StressMeasure::PK1: {os << "PK1"; break;} case StressMeasure::PK2: {os << "PK2"; break;} case StressMeasure::Kirchhoff: {os << "Kirchhoff"; break;} case StressMeasure::Biot: {os << "Biot"; break;} case StressMeasure::Mandel: {os << "Mandel"; break;} default: throw std::runtime_error ("a stress measure must be missing"); break; } return os; } /* ---------------------------------------------------------------------- */ std::ostream & operator<<(std::ostream & os, StrainMeasure s) { switch (s) { case StrainMeasure::Gradient: {os << "Gradient"; break;} case StrainMeasure::Infinitesimal: {os << "Infinitesimal"; break;} case StrainMeasure::GreenLagrange: {os << "Green-Lagrange"; break;} case StrainMeasure::Biot: {os << "Biot"; break;} case StrainMeasure::Log: {os << "Logarithmic"; break;} case StrainMeasure::Almansi: {os << "Almansi"; break;} case StrainMeasure::RCauchyGreen: {os << "Right Cauchy-Green"; break;} case StrainMeasure::LCauchyGreen: {os << "Left Cauchy-Green"; break;} default: throw std::runtime_error ("a strain measure must be missing"); } return os; } + /* ---------------------------------------------------------------------- */ + void banner(std::string name, Uint year, std::string cpy_holder) { + std::cout + << std::endl + << "µSpectre "<< name << std::endl + << "Copyright © " << year << " " << cpy_holder + << std::endl + << "This program comes with ABSOLUTELY NO WARRANTY." + << std::endl + << "This is free software, and you are welcome to redistribute it" + << std::endl + << "under certain conditions, see the license file." + << std::endl << std::endl; + } + } // muSpectre diff --git a/src/common/common.hh b/src/common/common.hh index acd0541..da04de4 100644 --- a/src/common/common.hh +++ b/src/common/common.hh @@ -1,201 +1,209 @@ /** * file common.hh * * @author Till Junge * * @date 01 May 2017 * * @brief Small definitions of commonly used types througout µSpectre * * @section LICENCE * * 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 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 #include #include #include #include +#include #ifndef COMMON_H #define COMMON_H namespace muSpectre { //! Eigen uses signed integers for dimensions. For consistency, µSpectre uses //! them througout the code using Dim_t = int;// needs to represent -1 for eigen constexpr Dim_t oneD{1}; constexpr Dim_t twoD{2}; constexpr Dim_t threeD{3}; constexpr Dim_t firstOrder{1}; constexpr Dim_t secondOrder{2}; constexpr Dim_t fourthOrder{4}; //! Scalar types used for mathematical calculations using Uint = unsigned int; using Int = int; using Real = double; using Complex = std::complex; //! Ccoord_t are cell coordinates, i.e. integer coordinates template using Ccoord_t = std::array; //! Real space coordinates template using Rcoord_t = std::array; template std::ostream & operator << (std::ostream & os, const std::array & index) { os << "("; for (size_t i = 0; i < dim-1; ++i) { os << index[i] << ", "; } os << index.back() << ")"; return os; } template Rcoord_t operator/(const Rcoord_t & a, const Rcoord_t & b) { Rcoord_t retval{a}; for (size_t i = 0; i < dim; ++i) { retval[i]/=b[i]; } return retval; } template Rcoord_t operator/(const Rcoord_t & a, const Ccoord_t & b) { Rcoord_t retval{a}; for (size_t i = 0; i < dim; ++i) { retval[i]/=b[i]; } return retval; } //! convenience definitions constexpr Real pi{3.1415926535897932384626433}; //! compile-time potentiation required for field-size computations template constexpr R ipow(R base, I exponent) { static_assert(std::is_integral::value, "Type must be integer"); R retval{1}; for (I i = 0; i < exponent; ++i) { retval *= base; } return retval; } + /** + * Copyright banner to be printed to the terminal by executables + * Arguments are the executable's name, year of writing and the name + * + address of the copyright holder + */ + void banner(std::string name, Uint year, std::string cpy_holder); + /** * Planner flags for FFT (follows FFTW, hopefully this choice will * be compatible with alternative FFT implementations) */ enum class FFT_PlanFlags {estimate, measure, patient}; //! continuum mechanics flags enum class Formulation{finite_strain, small_strain}; std::ostream & operator<<(std::ostream & os, Formulation f); /* ---------------------------------------------------------------------- */ //! Material laws can declare which type of stress measure they provide, //! and µSpectre will handle conversions enum class StressMeasure { Cauchy, PK1, PK2, Kirchhoff, Biot, Mandel, no_stress_}; std::ostream & operator<<(std::ostream & os, StressMeasure s); /* ---------------------------------------------------------------------- */ //! Material laws can declare which type of strain measure they require and //! µSpectre will provide it enum class StrainMeasure { Gradient, Infinitesimal, GreenLagrange, Biot, Log, Almansi, RCauchyGreen, LCauchyGreen, no_strain_}; std::ostream & operator<<(std::ostream & os, StrainMeasure s); /* ---------------------------------------------------------------------- */ /** Compile-time functions to set the stress and strain measures stored by mu_spectre depending on the formulation **/ constexpr StrainMeasure get_stored_strain_type(Formulation form) { switch (form) { case Formulation::finite_strain: { return StrainMeasure::Gradient; break; } case Formulation::small_strain: { return StrainMeasure::Infinitesimal; break; } default: return StrainMeasure::no_strain_; break; } } constexpr StressMeasure get_stored_stress_type(Formulation form) { switch (form) { case Formulation::finite_strain: { return StressMeasure::PK1; break; } case Formulation::small_strain: { return StressMeasure::Cauchy; break; } default: return StressMeasure::no_stress_; break; } } /* ---------------------------------------------------------------------- */ /** Compile-time functions to get the stress and strain measures after they may have been modified by choosing a formulation. For instance, a law that expecs a Green-Lagrange strain as input will get the infinitesimal strain tensor instead in a small strain computation **/ constexpr StrainMeasure get_formulation_strain_type(Formulation form, StrainMeasure expected) { switch (form) { case Formulation::finite_strain: { return expected; break; } case Formulation::small_strain: { return get_stored_strain_type(form); break; } default: return StrainMeasure::no_strain_; break; } } } // muSpectre #ifndef EXPLICITLY_TURNED_ON_CXX17 #include "common/utilities.hh" #endif #endif /* COMMON_H */ diff --git a/src/language_bindings/python/bind_py_common.cc b/src/language_bindings/python/bind_py_common.cc index b8d8fde..343f4af 100644 --- a/src/language_bindings/python/bind_py_common.cc +++ b/src/language_bindings/python/bind_py_common.cc @@ -1,115 +1,117 @@ /** * file bind_py_common.cc * * @author Till Junge * * @date 08 Jan 2018 * * @brief Python bindings for the common part of µSpectre * * @section LICENCE * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "common/ccoord_operations.hh" #include #include #include using namespace muSpectre; namespace py = pybind11; using namespace pybind11::literals; template void add_get_cube_helper(py::module & mod) { std::stringstream name {}; name << "get_" << dim << "d_cube"; mod.def (name.str().c_str(), &CcoordOps::get_cube, "size"_a, "return a Ccoord with the value 'size' repeated in each dimension"); } 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); } template void add_Pixels_helper(py::module & mod) { std::stringstream name{}; name << "Pixels" << dim << "d"; using Ccoord = 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); } PYBIND11_PLUGIN(common) { py::module mod("common", "Common utilities for pymuSpectre"); py::enum_(mod, "Formulation") .value("finite_strain", Formulation::finite_strain) //"µSpectre handles a problem in terms of tranformation gradient F and" //" first Piola-Kirchhoff stress P") .value("small_strain", Formulation::small_strain) //"µSpectre handles a problem in terms of the infinitesimal strain " //"tensor ε and Cauchy stress σ"); ; 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", FFT_PlanFlags::estimate) .value("measure", FFT_PlanFlags::measure) .value("patient", FFT_PlanFlags::patient); + mod.def("banner", &banner, "name"_a, "year"_a, "copyright_holder"_a); + add_get_cube(mod); add_Pixels(mod); return mod.ptr(); } diff --git a/src/language_bindings/python/bind_py_material.cc b/src/language_bindings/python/bind_py_material.cc index 5bfb052..db536b1 100644 --- a/src/language_bindings/python/bind_py_material.cc +++ b/src/language_bindings/python/bind_py_material.cc @@ -1,75 +1,89 @@ /** * file bind_py_material.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief python bindings for µSpectre's materials * * @section LICENCE * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "materials/material_hyper_elastic1.hh" +#include "system/system_base.hh" #include #include #include #include using namespace muSpectre; namespace py = pybind11; using namespace pybind11::literals; /** * python binding for the optionally objective form of Hooke's law */ template void add_material_hyper_elastic_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialHooke" << dim << 'd'; auto && name {name_stream.str().c_str()}; - py::class_>(mod, name) - .def(py::init(), "name"_a, "Young"_a, "Poisson"_a); + using Mat_t = MaterialHyperElastic1; + using Sys_t = SystemBase; + py::class_(mod, name) + .def(py::init(), "name"_a, "Young"_a, "Poisson"_a) + .def_static("make", + [](Sys_t & sys, std::string n, Real e, Real p) -> Mat_t & { + return Mat_t::make(sys, n, e, p); + }, + "system"_a, "name"_a, "Young"_a, "Poisson"_a, + py::return_value_policy::reference, py::keep_alive<1, 0>()) + .def("add_pixel", + [] (Mat_t & mat, Ccoord_t pix) { + mat.add_pixel(pix);}, + "pixel"_a); } + template void add_material_helper(py::module & mod) { add_material_hyper_elastic_helper(mod); } void add_material(py::module & mod) { add_material_helper(mod); add_material_helper(mod); } PYBIND11_PLUGIN(material) { (py::object) py::module::import("common"); (py::object) py::module::import("system"); py::module mod("material", "bindings for constitutive laws"); add_material(mod); return mod.ptr(); } diff --git a/src/language_bindings/python/bind_py_solvers.cc b/src/language_bindings/python/bind_py_solvers.cc index 28b1d1b..16b2f8b 100644 --- a/src/language_bindings/python/bind_py_solvers.cc +++ b/src/language_bindings/python/bind_py_solvers.cc @@ -1,139 +1,137 @@ /** * file bind_py_solver.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief python bindings for the muSpectre solvers * * @section LICENCE * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "solver/solvers.hh" #include #include #include #include using namespace muSpectre; namespace py=pybind11; using namespace pybind11::literals; /** * Solvers instanciated for systems with equal spatial and material dimension */ template void add_newton_cg_helper(py::module & mod) { std::stringstream name_stream {}; name_stream << "newton_cg" << sdim << "d"; - auto && name{name_stream.str().c_str()}; constexpr Dim_t mdim{sdim}; using sys = SystemBase; using grad = Grad_t; using grad_vec = Grad_t; - mod.def(name, + mod.def(name_stream.str().c_str(), [](sys & s, const grad & g, Formulation f, Real ct, Real nt, Uint max, Dim_t verb) -> typename sys::StrainField_t & { return newton_cg(s, g, f, ct, nt, max, verb); }, "system"_a, "ΔF₀"_a, "formulation"_a, "cg_tol"_a, "newton_tol"_a, "maxiter"_a=0, "verbose"_a=0); - mod.def(name, + mod.def(name_stream.str().c_str(), [](sys & s, const grad_vec & g, Formulation f, Real ct, Real nt, Uint max, Dim_t verb) -> typename sys::StrainField_t & { return newton_cg(s, g, f, ct, nt, max, verb); }, "system"_a, "ΔF₀"_a, "formulation"_a, "cg_tol"_a, "newton_tol"_a, "maxiter"_a=0, "verbose"_a=0); } template void add_de_geus_helper(py::module & mod) { std::stringstream name_stream {}; name_stream << "de_geus" << sdim << "d"; - auto && name{name_stream.str().c_str()}; constexpr Dim_t mdim{sdim}; using sys = SystemBase; using grad = Grad_t; using grad_vec = Grad_t; - mod.def(name, + mod.def(name_stream.str().c_str(), [](sys & s, const grad & g, Formulation f, Real ct, Real nt, Uint max, Dim_t verb) -> typename sys::StrainField_t & { return de_geus(s, g, f, ct, nt, max, verb); }, "system"_a, "ΔF₀"_a, "formulation"_a, "cg_tol"_a, "newton_tol"_a, "maxiter"_a=0, "verbose"_a=0); - mod.def(name, + mod.def(name_stream.str().c_str(), [](sys & s, const grad_vec & g, Formulation f, Real ct, Real nt, Uint max, Dim_t verb) -> typename sys::StrainField_t & { return de_geus(s, g, f, ct, nt, max, verb); }, "system"_a, "ΔF₀"_a, "formulation"_a, "cg_tol"_a, "newton_tol"_a, "maxiter"_a=0, "verbose"_a=0); } template void add_solver_helper(py::module & mod) { add_newton_cg_helper(mod); add_de_geus_helper (mod); } void add_solvers(py::module & mod) { add_solver_helper(mod); add_solver_helper(mod); } PYBIND11_PLUGIN(solvers) { py::module::import("common"); py::module::import("system"); py::module mod("solvers", "bindings for solvers"); add_solvers(mod); return mod.ptr(); } diff --git a/src/language_bindings/python/bind_py_system.cc b/src/language_bindings/python/bind_py_system.cc index 26ac23e..96f683a 100644 --- a/src/language_bindings/python/bind_py_system.cc +++ b/src/language_bindings/python/bind_py_system.cc @@ -1,95 +1,100 @@ /** * file bind_py_system.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief Python bindings for the system factory function * * @section LICENCE * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "common/ccoord_operations.hh" #include "system/system_factory.hh" #include "system/system_base.hh" #include #include #include using namespace muSpectre; namespace py=pybind11; using namespace pybind11::literals; /** * the system factory is only bound for default template parameters */ template void add_system_factory_helper(py::module & mod) { std::stringstream name {}; name << "SystemFactory" << dim << "d"; using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; using Form = Formulation; mod.def (name.str().c_str(), [](Ccoord res, Rcoord lens, Form form) { return make_system(std::move(res), std::move(lens), std::move(form)); }, "resolutions"_a, "lengths"_a=CcoordOps::get_cube(1.), "formulation"_a=Formulation::finite_strain); } void add_system_factory(py::module & mod) { - add_system_factory_helper(mod); + add_system_factory_helper(mod); add_system_factory_helper(mod); } /** * SystemBase for which the material and spatial dimension are identical */ template void add_system_base_helper(py::module & mod) { std::stringstream name{}; name << "SystemBase" << dim << 'd'; auto && name_str{name.str().c_str()}; - - py::class_>(mod, name_str); + using sys_t = SystemBase; + py::class_(mod, name_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=FFT_PlanFlags::estimate); } void add_system_base(py::module & mod) { add_system_base_helper (mod); add_system_base_helper(mod); } PYBIND11_PLUGIN(system) { py::module::import("common"); py::module mod("system", "bindings for systems and system factories"); add_system_factory(mod); add_system_base(mod); return mod.ptr(); } diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index f9d8d73..0e0c546 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,637 +1,661 @@ /** * file material_muSpectre_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials written for µSpectre specifically. These * can take full advantage of the configuration-change utilities of * µSpectre. The user can inherit from them to define new constitutive * laws and is merely required to provide the methods for computing the * second Piola-Kirchhoff stress and Tangent. This class uses the * "curiously recurring template parameter" to avoid virtual calls. * * @section LICENCE * * 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 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. */ #ifndef MATERIAL_MUSPECTRE_BASE_H #define MATERIAL_MUSPECTRE_BASE_H #include "common/common.hh" #include "materials/material_base.hh" #include "materials/materials_toolbox.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "common//utilities.hh" #include #include #include #include namespace muSpectre { + /** + * Forward declaration for factory function + */ + template + class SystemBase; + template struct MaterialMuSpectre_traits { }; template class MaterialMuSpectre; template <> struct MaterialMuSpectre_traits { using DefaultInternalVariables = std::tuple<>; }; //! 'Material' is a CRTP template class MaterialMuSpectre: public MaterialBase { public: using NeedTangent = MatTB::NeedTangent; using Parent = MaterialBase; using GFieldCollection_t = typename Parent::GFieldCollection_t; using MFieldCollection_t = typename Parent::MFieldCollection_t; using StressField_t = typename Parent::StressField_t; using StrainField_t = typename Parent::StrainField_t; using TangentField_t = typename Parent::TangentField_t; using DefaultInternalVariables = std::tuple<>; using traits = MaterialMuSpectre_traits; //! Default constructor MaterialMuSpectre() = delete; //! Construct by name MaterialMuSpectre(std::string name); //! Copy constructor MaterialMuSpectre(const MaterialMuSpectre &other) = delete; //! Move constructor MaterialMuSpectre(MaterialMuSpectre &&other) noexcept = delete; //! Destructor virtual ~MaterialMuSpectre() noexcept = default; + //! Factory + template + static Material & make(SystemBase & sys, + ConstructorArgs &&... args); + //! Copy assignment operator MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete; //! Move assignment operator MaterialMuSpectre& operator=(MaterialMuSpectre &&other) noexcept = delete; //* allocate memory, etc virtual void initialise(bool stiffness = false) override final; using Parent::compute_stresses; using Parent::compute_stresses_tangent; //! computes stress virtual void compute_stresses(const StrainField_t & F, StressField_t & P, Formulation form) override final; //! computes stress and tangent modulus virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form) override final; protected: //! computes stress with the formulation available at compile time //! __attribute__ required by g++-6 and g++-7 because of this bug: //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947 template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P) __attribute__ ((visibility ("default"))); //! computes stress with the formulation available at compile time //! __attribute__ required by g++-6 and g++-7 because of this bug: //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947 template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P, TangentField_t & K) __attribute__ ((visibility ("default"))); //! this iterable class is a default for simple laws that just take a strain //! the iterable is just a templated wrapper to provide a range to iterate over //! that does or does not include tangent moduli template class iterable_proxy; /** * inheriting classes with internal variables need to overload this function */ decltype(auto) get_internals() { // the default material has no internal variables return typename Material::InternalVariables{};} decltype(auto) get_internals() const { // the default material has no internal variables return typename Material::InternalVariables{};} typename traits::InternalVariables internal_variables{}; bool is_initialised{false}; private: }; /* ---------------------------------------------------------------------- */ template MaterialMuSpectre:: MaterialMuSpectre(std::string name) :Parent(name) { using stress_compatible = typename Material::StressMap_t:: template is_compatible; using strain_compatible = typename Material::StrainMap_t:: template is_compatible; using tangent_compatible = typename Material::TangentMap_t:: template is_compatible; static_assert((stress_compatible::value && stress_compatible::explain()), "The material's declared stress map is not compatible " "with the stress field. More info in previously shown " "assert."); static_assert((strain_compatible::value && strain_compatible::explain()), "The material's declared strain map is not compatible " "with the strain field. More info in previously shown " "assert."); static_assert((tangent_compatible::value && tangent_compatible::explain()), "The material's declared tangent map is not compatible " "with the tangent field. More info in previously shown " "assert."); } + /* ---------------------------------------------------------------------- */ + template + template + Material & MaterialMuSpectre:: + make(SystemBase & sys, + ConstructorArgs && ... args) { + auto mat = std::make_unique(args...); + auto & mat_ref = *mat; + sys.add_material(std::move(mat)); + return mat_ref; + } + + /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: initialise(bool /*stiffness*/) { if (!this->is_initialised) { this->internal_fields.initialise(); this->is_initialised = true; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: compute_stresses(const StrainField_t &F, StressField_t &P, Formulation form) { switch (form) { case Formulation::finite_strain: { this->template compute_stresses_worker(F, P); break; } case Formulation::small_strain: { this->template compute_stresses_worker(F, P); break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form) { switch (form) { case Formulation::finite_strain: { this->template compute_stresses_worker(F, P, K); break; } case Formulation::small_strain: { this->template compute_stresses_worker(F, P, K); break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre:: compute_stresses_worker(const StrainField_t & F, StressField_t & P, TangentField_t & K){ /* These lambdas are executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ using Strains_t = std::tuple; using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this] (Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, Material::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto && F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // return value contains a tuple of rvalue_refs to both stress and tangent moduli Stresses = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); }; auto constitutive_law_finite_strain = [this] (Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, Material::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto && F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // TODO: Figure this out: I can't std::move(internals...), // because if there are no internals, compilation fails with "no // matching function for call to ‘move()’'. These are tuples of // lvalue references, so it shouldn't be too bad, but still // irksome. // return value contains a tuple of rvalue_refs to both stress // and tangent moduli auto stress_tgt = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); auto && stress = std::get<0>(stress_tgt); auto && tangent = std::get<1>(stress_tgt); Stresses = MatTB::PK1_stress (std::move(F), std::move(stress), std::move(tangent)); }; iterable_proxy fields{*this, F, P, K}; for (auto && arglist: fields) { /** * arglist is a tuple of three tuples containing only Lvalue * references (see value_tye in the class definition of * iterable_proxy::iterator). Tuples contain strains, stresses * and internal variables, respectively, */ //auto && stress_tgt = std::get<0>(tuples); //auto && inputs = std::get<1>(tuples);TODO:clean this static_assert(std::is_same(std::get<0>(arglist)))>>::value, "Type mismatch for strain reference, check iterator " "value_type"); static_assert(std::is_same(std::get<1>(arglist)))>>::value, "Type mismatch for stress reference, check iterator" "value_type"); static_assert(std::is_same(std::get<1>(arglist)))>>::value, "Type mismatch for tangent reference, check iterator" "value_type"); switch (Form) { case Formulation::small_strain: { apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre:: compute_stresses_worker(const StrainField_t & F, StressField_t & P){ /* These lambdas are executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ using Strains_t = std::tuple; using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this] (Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, Material::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto && F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // return value contains a tuple of rvalue_refs to both stress and tangent moduli auto && sigma = std::get<0>(Stresses); sigma = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress(std::move(strain), internals...);}, internal_variables); }; auto constitutive_law_finite_strain = [this] (Strains_t Strains, Stresses_t && Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, Material::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto && F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // TODO: Figure this out: I can't std::move(internals...), // because if there are no internals, compilation fails with "no // matching function for call to ‘move()’'. These are tuples of // lvalue references, so it shouldn't be too bad, but still // irksome. // return value contains a tuple of rvalue_refs to both stress // and tangent moduli auto && stress = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress(std::move(strain), internals...);}, internal_variables); auto && P = get<0>(Stresses); P = MatTB::PK1_stress (F, stress); }; iterable_proxy fields{*this, F, P}; for (auto && arglist: fields) { /** * arglist is a tuple of three tuples containing only Lvalue * references (see value_tye in the class definition of * iterable_proxy::iterator). Tuples contain strains, stresses * and internal variables, respectively, */ //auto && stress_tgt = std::get<0>(tuples); //auto && inputs = std::get<1>(tuples);TODO:clean this static_assert(std::is_same(std::get<0>(arglist)))>>::value, "Type mismatch for strain reference, check iterator " "value_type"); static_assert(std::is_same(std::get<1>(arglist)))>>::value, "Type mismatch for stress reference, check iterator" "value_type"); switch (Form) { case Formulation::small_strain: { apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ //! this iterator class is a default for simple laws that just take a strain template template class MaterialMuSpectre::iterable_proxy { public: //! Default constructor iterable_proxy() = delete; using NeedTangent = typename MaterialMuSpectre::NeedTangent; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ template iterable_proxy(const MaterialMuSpectre & mat, const StrainField_t & F, StressField_t & P, std::enable_if_t & K) :material(mat), strain_field(F), stress_tup{P,K}, internals(material.internal_variables){}; template iterable_proxy(const MaterialMuSpectre & mat, const StrainField_t & F, std::enable_if_t & P) :material(mat), strain_field(F), stress_tup{P}, internals(material.internal_variables){}; using StrainMap_t = typename Material::StrainMap_t; using StressMap_t = typename Material::StressMap_t; using TangentMap_t = typename Material::TangentMap_t; using Strain_t = typename Material::Strain_t; using Stress_t = typename Material::Stress_t; using Tangent_t = typename Material::Tangent_t; using InternalVariables = typename Material::InternalVariables; using StressFieldTup = std::conditional_t <(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; using StressMapTup = std::conditional_t <(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; using Stress_tTup = std::conditional_t<(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! Copy constructor iterable_proxy(const iterable_proxy &other) = default; //! Move constructor iterable_proxy(iterable_proxy &&other) noexcept = default; //! Destructor virtual ~iterable_proxy() noexcept = default; //! Copy assignment operator iterable_proxy& operator=(const iterable_proxy &other) = default; //! Move assignment operator iterable_proxy& operator=(iterable_proxy &&other) = default; class iterator { public: using InternalReferences = MatTB::ReferenceTuple_t; using value_type = std::tuple, Stress_tTup, InternalReferences>; using iterator_category = std::forward_iterator_tag; //! Default constructor iterator() = delete; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ iterator(const iterable_proxy & it, bool begin = true) : it{it}, strain_map{it.strain_field}, stress_map {it.stress_tup}, index{begin ? 0:it.material.internal_fields.size()}{} //! Copy constructor iterator(const iterator &other) = default; //! Move constructor iterator(iterator &&other) noexcept = default; //! Destructor virtual ~iterator() noexcept = default; //! Copy assignment operator iterator& operator=(const iterator &other) = default; //! Move assignment operator iterator& operator=(iterator &&other) = default; //! pre-increment inline iterator & operator++(); //! dereference inline value_type operator*(); //! inequality inline bool operator!=(const iterator & other) const; protected: const iterable_proxy & it; StrainMap_t strain_map; StressMapTup stress_map; size_t index; private: }; iterator begin() {return std::move(iterator(*this));} iterator end() {return std::move(iterator(*this, false));} protected: const MaterialMuSpectre & material; const StrainField_t & strain_field; StressFieldTup stress_tup; const InternalVariables & internals; private: }; /* ---------------------------------------------------------------------- */ template template bool MaterialMuSpectre::iterable_proxy::iterator:: operator!=(const iterator & other) const { return (this->index != other.index); } /* ---------------------------------------------------------------------- */ template template typename MaterialMuSpectre:: template iterable_proxy:: iterator & MaterialMuSpectre::iterable_proxy::iterator:: operator++() { this->index++; return *this; } /* ---------------------------------------------------------------------- */ template template typename MaterialMuSpectre:: template iterable_proxy::iterator:: value_type MaterialMuSpectre::iterable_proxy::iterator:: operator*() { const Ccoord_t pixel{ this->it.material.internal_fields.get_ccoord(this->index)}; auto && strain = std::make_tuple(this->strain_map[pixel]); auto && stresses = apply([&pixel] (auto && ... stress_tgt) { return std::make_tuple(stress_tgt[pixel]...);}, this->stress_map); const auto & internal = this->it.material.get_internals(); const auto id{index}; auto && internals = apply([id] (auto && ... internals) { return std::make_tuple(internals[id]...);}, internal); return std::make_tuple(std::move(strain), std::move(stresses), std::move(internals)); } } // muSpectre #endif /* MATERIAL_MUSPECTRE_BASE_H */ diff --git a/src/system/system_base.cc b/src/system/system_base.cc index 1a8555e..de159f0 100644 --- a/src/system/system_base.cc +++ b/src/system/system_base.cc @@ -1,202 +1,204 @@ /** * file system_base.cc * * @author Till Junge * * @date 01 Nov 2017 * * @brief Implementation for system base class * * @section LICENCE * * 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 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 "system/system_base.hh" #include "common/ccoord_operations.hh" #include "common/iterators.hh" #include "common/tensor_algebra.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template SystemBase::SystemBase(Projection_ptr projection_, Formulation form) :resolutions{projection_->get_resolutions()}, pixels(resolutions), lengths{projection_->get_lengths()}, fields{}, F{make_field("Gradient", this->fields)}, P{make_field("Piola-Kirchhoff-1", this->fields)}, projection{std::move(projection_)}, form{form} { } /* ---------------------------------------------------------------------- */ template - void SystemBase::add_material(Material_ptr mat) { + typename SystemBase::Material_t & + SystemBase::add_material(Material_ptr mat) { this->materials.push_back(std::move(mat)); + return *this->materials.back(); } /* ---------------------------------------------------------------------- */ template typename SystemBase::FullResponse_t SystemBase::evaluate_stress_tangent(StrainField_t & grad) { if (this->is_initialised == false) { this->initialise(); } //! High level compatibility checks if (grad.size() != this->F.size()) { throw std::runtime_error("Size mismatch"); } if (!this->K) { this->K = make_field("Tangent Stiffness", this->fields); } for (auto & mat: this->materials) { mat->compute_stresses_tangent(grad, this->P, this->K.value(), this->form); } return std::tie(this->P, this->K.value()); } /* ---------------------------------------------------------------------- */ template typename SystemBase::StressField_t & SystemBase::directional_stiffness(const TangentField_t &K, const StrainField_t &delF, StressField_t &delP) { for (auto && tup: akantu::zip(K.get_map(), delF.get_map(), delP.get_map())){ auto & k = std::get<0>(tup); auto & df = std::get<1>(tup); auto & dp = std::get<2>(tup); dp = Matrices::tensmult(k, df); } return this->project(delP); } /* ---------------------------------------------------------------------- */ template typename SystemBase::StressField_t & SystemBase::project(StressField_t &field) { this->projection->apply_projection(field); return field; } /* ---------------------------------------------------------------------- */ template typename SystemBase::StrainField_t & SystemBase::get_strain() { if (this->is_initialised == false) { this->initialise(); } return this->F; } /* ---------------------------------------------------------------------- */ template const typename SystemBase::StressField_t & SystemBase::get_stress() const { return this->P; } /* ---------------------------------------------------------------------- */ template void SystemBase::initialise(FFT_PlanFlags flags) { // check that all pixels have been assigned exactly one material this->check_material_coverage(); // resize all global fields (strain, stress, etc) this->fields.initialise(this->resolutions); // initialise the projection and compute the fft plan this->projection->initialise(flags); this->is_initialised = true; } /* ---------------------------------------------------------------------- */ template void SystemBase::initialise_materials(bool stiffness) { for (auto && mat: this->materials) { mat->initialise(stiffness); } } /* ---------------------------------------------------------------------- */ template typename SystemBase::iterator SystemBase::begin() { return this->pixels.begin(); } /* ---------------------------------------------------------------------- */ template typename SystemBase::iterator SystemBase::end() { return this->pixels.end(); } /* ---------------------------------------------------------------------- */ template void SystemBase::check_material_coverage() { auto nb_pixels = CcoordOps::get_size(this->resolutions); std::vector*> assignments(nb_pixels, nullptr); for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->resolutions, pixel); auto& assignment{assignments.at(index)}; if (assignment != nullptr) { std::stringstream err{}; err << "Pixel " << pixel << "is already assigned to material '" << assignment->get_name() << "' and cannot be reassigned to material '" << mat->get_name(); throw std::runtime_error(err.str()); } else { assignments[index] = mat.get(); } } } // find and identify unassigned pixels std::vector unassigned_pixels; for (size_t i = 0; i < assignments.size(); ++i) { if (assignments[i] == nullptr) { unassigned_pixels.push_back(CcoordOps::get_ccoord(this->resolutions, i)); } } if (unassigned_pixels.size() != 0) { std::stringstream err {}; err << "The following pixels have were not assigned a material: "; for (auto & pixel: unassigned_pixels) { err << pixel << ", "; } err << "and that cannot be handled"; throw std::runtime_error(err.str()); } } template class SystemBase; template class SystemBase; } // muSpectre diff --git a/src/system/system_base.hh b/src/system/system_base.hh index c4f09ff..fa5bd1e 100644 --- a/src/system/system_base.hh +++ b/src/system/system_base.hh @@ -1,160 +1,160 @@ /** * file system_base.hh * * @author Till Junge * * @date 01 Nov 2017 * * @brief Base class representing a unit cell system with single * projection operator * * @section LICENCE * * 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 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. */ #ifndef SYSTEM_BASE_H #define SYSTEM_BASE_H #include "common/common.hh" #include "common/ccoord_operations.hh" #include "common/field.hh" #include "common/utilities.hh" #include "materials/material_base.hh" #include "fft/projection_base.hh" #include #include #include #include namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) - template + template class SystemBase { public: using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; using FieldCollection_t = FieldCollection; using Material_t = MaterialBase; using Material_ptr = std::unique_ptr; using Projection_ptr = std::unique_ptr>; using StrainField_t = TensorField; using StressField_t = TensorField; using TangentField_t = TensorField; using FullResponse_t = std::tuple; using iterator = typename CcoordOps::Pixels::iterator; //! Default constructor SystemBase() = delete; //! constructor using sizes and resolution SystemBase(Projection_ptr projection, Formulation form=Formulation::finite_strain); //! Copy constructor SystemBase(const SystemBase &other) = delete; //! Move constructor SystemBase(SystemBase &&other) noexcept = default; //! Destructor virtual ~SystemBase() noexcept = default; //! Copy assignment operator SystemBase& operator=(const SystemBase &other) = delete; //! Move assignment operator SystemBase& operator=(SystemBase &&other) = default; /** * Materials can only be moved. This is to assure exclusive * ownership of any material by this system */ - void add_material(Material_ptr mat); + Material_t & add_material(Material_ptr mat); /** * evaluates all materials */ FullResponse_t evaluate_stress_tangent(StrainField_t & F); /** * evaluate directional stiffness (i.e. G:K:δF) */ StressField_t & directional_stiffness(const TangentField_t & K, const StrainField_t & delF, StressField_t & delP); /** * Convenience function circumventing the neeed to use the * underlying projection */ StressField_t & project(StressField_t & field); StrainField_t & get_strain(); const StressField_t & get_stress() const; /** * general initialisation; initialises the projection and * fft_engine (i.e. infrastructure) but not the materials. These * need to be initialised separately */ void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate); /** * initialise materials (including resetting any history variables) */ void initialise_materials(bool stiffness=false); iterator begin(); iterator end(); size_t size() const {return pixels.size();} const Ccoord & get_resolutions() const {return this->resolutions;} const Rcoord & get_lengths() const {return this->lengths;} protected: //! make sure that every pixel is assigned to one and only one material void check_material_coverage(); const Ccoord & resolutions; CcoordOps::Pixels pixels; const Rcoord & lengths; FieldCollection_t fields; StrainField_t & F; StressField_t & P; //! Tangent field might not even be required; so this is an //! optional ref_wrapper instead of a ref optional> K{}; std::vector materials{}; Projection_ptr projection; Formulation form; bool is_initialised{false}; private: }; } // muSpectre #endif /* SYSTEM_BASE_H */ diff --git a/tests/python_binding_tests.py b/tests/python_binding_tests.py index b97d84d..99f1b99 100755 --- a/tests/python_binding_tests.py +++ b/tests/python_binding_tests.py @@ -1,60 +1,113 @@ #!/usr/bin/env python3 """ file python_binding_tests.py @author Till Junge @date 09 Jan 2018 @brief Unit tests for python bindings @section LICENCE 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. """ import unittest import sys +import os +import numpy as np -sys.path.append("src/language_bindings/python") +sys.path.append(os.path.join(os.getcwd(), "src/language_bindings/python")) try: import common import system import solvers import material except ImportError as err: print(err) sys.exit(-1) class SystemCheck(unittest.TestCase): def test_Construction(self): """ Simple check for system constructors """ - print(dir(system)) - print(help(system.SystemFactory2d)) + resolution = [5,7] + lengths = [5.2, 8.3] + formulation = common.Formulation.small_strain + try: + sys = system.SystemFactory2d(resolution, + lengths, + formulation) + mat = material.MaterialHooke2d.make(sys, "material", 210e9, .33) + except Exception(err): + print(err) + raise err + +class MaterialHooke2dCheck(unittest.TestCase): + def setUp(self): + self.resolution = [5,7] + self.lengths = [5.2, 8.3] + self.formulation = common.Formulation.small_strain + self.sys = system.SystemFactory2d(self.resolution, + self.lengths, + self.formulation) + self.mat = material.MaterialHooke2d.make( + self.sys, "material", 210e9, .33) + + def test_add_material(self): + self.mat.add_pixel([2,1]) + class SolverCheck(unittest.TestCase): - pass + def setUp(self): + self.resolution = [5,7] + self.lengths = [5.2, 8.3] + self.formulation = common.Formulation.small_strain + self.sys = system.SystemFactory2d(self.resolution, + self.lengths, + self.formulation) + self.hard = material.MaterialHooke2d.make( + self.sys, "hard", 210e9, .33) + self.soft = material.MaterialHooke2d.make( + self.sys, "soft", 70e9, .33) + + def test_solve(self): + for i, pixel in enumerate(self.sys): + if i < 3: + self.hard.add_pixel(pixel) + else: + self.soft.add_pixel(pixel) + + self.sys.initialise() + tol = 1e-6 + Del0 = np.array([[0, .1], + [0, 0]]) + maxiter = 100 + verbose = 2 + # the following segfaults: + #TODO solvers.de_geus2d(self.sys, Del0, self.formulation, + # tol, tol, maxiter, verbose) if __name__ == '__main__': unittest.main() diff --git a/tests/test_solver_newton_cg.cc b/tests/test_solver_newton_cg.cc index fcfb7ed..54975d5 100644 --- a/tests/test_solver_newton_cg.cc +++ b/tests/test_solver_newton_cg.cc @@ -1,175 +1,174 @@ /** * file test_solver_newton_cg.cc * * @author Till Junge * * @date 20 Dec 2017 * * @brief Tests for the standard Newton-Raphson + Conjugate Gradient solver * * @section LICENCE * * 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 GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "tests.hh" #include "solver/solvers.hh" #include "fft/fftw_engine.hh" #include "fft/projection_finite_strain_fast.hh" #include "materials/material_hyper_elastic1.hh" #include "common/iterators.hh" #include "common/ccoord_operations.hh" #include "system/system_factory.hh" 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, lengths)}; auto proj_ptr{std::make_unique>(std::move(fft_ptr))}; SystemBase sys(std::move(proj_ptr)); using Mat_t = MaterialHyperElastic1; //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 = std::make_unique("hard", 10*Young, Poisson); - auto Material_soft = std::make_unique("soft", Young, 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); + Material_hard.add_pixel(pixel); } else { - Material_soft->add_pixel(pixel); + Material_soft.add_pixel(pixel); } } - sys.add_material(std::move(Material_hard)); - sys.add_material(std::move(Material_soft)); 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 bool verbose{false}; constexpr Formulation form{Formulation::finite_strain}; GradIncrements grads; grads.push_back(delF0); Eigen::ArrayXXd res1{de_geus(sys, grads, form, cg_tol, newton_tol, maxiter, verbose).eigen()}; Eigen::ArrayXXd res2{newton_cg(sys, grads, form, cg_tol, newton_tol, maxiter, verbose).eigen()}; 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{3, 3}; constexpr Rcoord lengths{3, 3}; 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_system(resolutions, lengths, form)}; using Mat_t = MaterialHyperElastic1; 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; constexpr Real eps0 = 1.; delEps0 << eps0, 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 bool verbose{false}; auto & result = newton_cg(sys, delEps0, form, cg_tol, newton_tol, maxiter, verbose); if (verbose) { std::cout << "result:" << std::endl << result.eigen() << std::endl; std::cout << "mean strain = " << std::endl << result.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; for (const auto & pixel: sys) { if (pixel[0] < Dim_t(nb_lays)) { BOOST_CHECK_LE((Eps_hard-result.get_map()[pixel]).norm(), tol); } else { BOOST_CHECK_LE((Eps_soft-result.get_map()[pixel]).norm(), tol); } } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre