diff --git a/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc index 0e09ef5..6af8805 100644 --- a/language_bindings/python/bind_py_material.cc +++ b/language_bindings/python/bind_py_material.cc @@ -1,252 +1,317 @@ /** * @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 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_linear_elastic1.hh" #include "materials/material_linear_elastic2.hh" #include "materials/material_linear_elastic3.hh" #include "materials/material_linear_elastic4.hh" +#include "materials/material_anisotropic.hh" +#include "materials/material_orthotropic.hh" #include "materials/material_base.hh" #include "cell/cell_base.hh" #include #include #include #include #include using namespace muSpectre; namespace py = pybind11; using namespace pybind11::literals; template void add_material_base_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialBase" << dim << 'd'; const auto name {name_stream.str()}; using Mat_t = MaterialBase; py::class_(mod, name.c_str()); } +/** + * python binding for the optionally objective form of Anisotropic material + */ +template +void add_material_anisotropic_helper(py::module & mod) { + std::stringstream name_stream{}; + name_stream << "MaterialAnisotropic_" << dim << 'd'; + const auto name {name_stream.str()}; + using Mat_t = MaterialAnisotropic; + using Sys_t = CellBase; + using MatBase_t = MaterialBase; + py::class_(mod, name.c_str()) + .def_static("make", + [](Sys_t & sys, std::string n, std::vector input) -> Mat_t & { + return Mat_t::make(sys, n, input); + }, + "cell"_a, "name"_a, "inputs"_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) + .def("add_pixel_split", + [] (Mat_t & mat, Ccoord_t pix, Real ratio) { + mat.add_pixel_split(pix, ratio);}, + "pixel"_a, + "ratio"_a) + .def("size", &Mat_t::size); +} + + +/** + * python binding for the optionally objective form of Anisotropic material + */ +template +void add_material_orthotropic_helper(py::module & mod) { + std::stringstream name_stream{}; + name_stream << "MaterialOrthotropic_" << dim << 'd'; + const auto name {name_stream.str()}; + using Mat_t = MaterialOrthotropic; + using Sys_t = CellBase; + using MatAniso_t = MaterialAnisotropic; + // using MatBase_t = MaterialBase; + py::class_(mod, name.c_str()) + .def_static("make", + [](Sys_t & sys, std::string n, std::vector input) -> Mat_t & { + return Mat_t::make(sys, n, input); + }, + "cell"_a, "name"_a, "inputs"_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) + .def("add_pixel_split", + [] (Mat_t & mat, Ccoord_t pix, Real ratio) { + mat.add_pixel_split(pix, ratio);}, + "pixel"_a, + "ratio"_a) + .def("size", &Mat_t::size); + } /** * python binding for the optionally objective form of Hooke's law */ template void add_material_linear_elastic1_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic1_" << dim << 'd'; const auto name {name_stream.str()}; using Mat_t = MaterialLinearElastic1; using Sys_t = CellBase; using MatBase_t = MaterialBase; py::class_(mod, name.c_str()) .def_static("make", [](Sys_t & sys, std::string n, Real e, Real p) -> Mat_t & { return Mat_t::make(sys, n, e, p); }, "cell"_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) .def("add_pixel_split", [] (Mat_t & mat, Ccoord_t pix, Real ratio) { mat.add_pixel_split(pix, ratio);}, "pixel"_a, "ratio"_a) .def("size", &Mat_t::size); } template void add_material_linear_elastic2_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic2_" << dim << 'd'; const auto name {name_stream.str()}; using Mat_t = MaterialLinearElastic2; using Sys_t = CellBase; using MatBase_t = MaterialBase; py::class_(mod, name.c_str()) .def_static("make", [](Sys_t & sys, std::string n, Real e, Real p) -> Mat_t & { return Mat_t::make(sys, n, e, p); }, "cell"_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, py::EigenDRef& eig) { Eigen::Matrix eig_strain{eig}; mat.add_pixel(pix, eig_strain);}, "pixel"_a, "eigenstrain"_a) .def("add_pixel_split", [] (Mat_t & mat, Ccoord_t pix, Real ratio, py::EigenDRef& eig) { Eigen::Matrix eig_strain{eig}; mat.add_pixel_split(pix, ratio, eig_strain);}, "pixel"_a, "ratio"_a, "eigenstrain"_a) .def("size", &Mat_t::size); } template void add_material_linear_elastic3_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic3_" << dim << 'd'; const auto name {name_stream.str()}; using Mat_t = MaterialLinearElastic3; using Sys_t = CellBase; using MatBase_t = MaterialBase; py::class_(mod, name.c_str()) .def(py::init(), "name"_a) .def_static("make", [](Sys_t & sys, std::string n) -> Mat_t & { return Mat_t::make(sys, n); }, "cell"_a, "name"_a, py::return_value_policy::reference, py::keep_alive<1, 0>()) .def("add_pixel", [] (Mat_t & mat, Ccoord_t pix, Real Young, Real Poisson) { mat.add_pixel(pix, Young, Poisson);}, "pixel"_a, "Young"_a, "Poisson"_a) .def("add_pixel_split", [] (Mat_t & mat, Ccoord_t pix, Real ratio, Real Young, Real Poisson) { mat.add_pixel_split(pix, ratio, Young, Poisson);}, "pixel"_a, "ratio"_a, "Young"_a, "Poisson"_a) .def("size", &Mat_t::size); } template void add_material_linear_elastic4_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic4_" << dim << 'd'; const auto name {name_stream.str()}; using Mat_t = MaterialLinearElastic4; using Sys_t = CellBase; py::class_(mod, name.c_str()) .def(py::init(), "name"_a) .def_static("make", [](Sys_t & sys, std::string n) -> Mat_t & { return Mat_t::make(sys, n); }, "cell"_a, "name"_a, py::return_value_policy::reference, py::keep_alive<1, 0>()) .def("add_pixel", [] (Mat_t & mat, Ccoord_t pix, Real Young, Real Poisson) { mat.add_pixel(pix, Young, Poisson);}, "pixel"_a, "Young"_a, "Poisson"_a) .def("add_pixel_split", [] (Mat_t & mat, Ccoord_t pix, Real ratio, Real Young, Real Poisson) { mat.add_pixel_split(pix, ratio, Young, Poisson);}, "pixel"_a, "ratio"_a, "Young"_a, "Poisson"_a) .def("size", &Mat_t::size); } template class PyMaterialBase : public MaterialBase { public: /* Inherit the constructors */ using Parent = 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) */ ); } virtual void compute_stresses(const typename Parent::StrainField_t & F, typename Parent::StressField_t & P, 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 ); } virtual void compute_stresses_tangent(const typename Parent::StrainField_t & F, typename Parent::StressField_t & P, typename Parent::TangentField_t & K, Formulation form) override { PYBIND11_OVERLOAD_PURE (void, /* Return type */ Parent, /* Parent class */ compute_stresses_tangent, /* Name of function in C++ (must match Python name) */ F,P,K, form ); } }; template void add_material_helper(py::module & mod) { add_material_base_helper(mod); + add_material_anisotropic_helper(mod); + add_material_orthotropic_helper(mod); add_material_linear_elastic1_helper(mod); add_material_linear_elastic2_helper(mod); add_material_linear_elastic3_helper(mod); add_material_linear_elastic4_helper(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); } diff --git a/language_bindings/python/bind_py_material.cc.rej b/language_bindings/python/bind_py_material.cc.rej new file mode 100644 index 0000000..a9cd953 --- /dev/null +++ b/language_bindings/python/bind_py_material.cc.rej @@ -0,0 +1,30 @@ +diff a/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc (rejected hunks) +@@ -46,15 +47,26 @@ using namespace pybind11::literals; + /** + * python binding for the optionally objective form of Hooke's law + */ ++template ++void add_material_base_helper(py::module & mod) { ++ std::stringstream name_stream{}; ++ name_stream << "MaterialBase_" << dim << 'd'; ++ const auto name {name_stream.str()}; ++ ++ using Mat_t = MaterialBase; ++ //using Sys_t = CellBase; ++ py::class_(mod, name.c_str()); ++} ++ + template + void add_material_linear_elastic1_helper(py::module & mod) { + std::stringstream name_stream{}; + name_stream << "MaterialLinearElastic1_" << dim << 'd'; + const auto name {name_stream.str()}; +- ++ using MatBase_t = MaterialBase; + using Mat_t = MaterialLinearElastic1; + using Sys_t = CellBase; +- py::class_(mod, name.c_str()) ++ py::class_(mod, name.c_str()) + .def_static("make", + [](Sys_t & sys, std::string n, Real e, Real p) -> Mat_t & { + return Mat_t::make(sys, n, e, p); diff --git a/src/cell/cell_base.cc b/src/cell/cell_base.cc index 9f4f7bf..a5578dd 100644 --- a/src/cell/cell_base.cc +++ b/src/cell/cell_base.cc @@ -1,485 +1,484 @@ /* * @file cell_base.cc * * @author Till Junge * * @date 01 Nov 2017 * * @brief Implementation for cell base class * * 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 "cell/cell_base.hh" #include "common/ccoord_operations.hh" #include "common/iterators.hh" #include "common/tensor_algebra.hh" #include "common/common.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template - CellBase::CellBase(Projection_ptr projection_, SplittedCell is_cell_splitted) + CellBase::CellBase(Projection_ptr projection_) :subdomain_resolutions{projection_->get_subdomain_resolutions()}, subdomain_locations{projection_->get_subdomain_locations()}, domain_resolutions{projection_->get_domain_resolutions()}, pixels(subdomain_resolutions, subdomain_locations), domain_lengths{projection_->get_domain_lengths()}, fields{std::make_unique()}, F{make_field("Gradient", *this->fields)}, P{make_field("Piola-Kirchhoff-1", *this->fields)}, - projection{std::move(projection_)}, - is_cell_splitted{is_cell_splitted} + projection{std::move(projection_)} { } /** * turns out that the default move container in combination with * clang segfaults under certain (unclear) cicumstances, because the * move constructor of the optional appears to be busted in gcc * 7.2. Copying it (K) instead of moving it fixes the issue, and * since it is a reference, the cost is practically nil */ template CellBase::CellBase(CellBase && other): subdomain_resolutions{std::move(other.subdomain_resolutions)}, subdomain_locations{std::move(other.subdomain_locations)}, domain_resolutions{std::move(other.domain_resolutions)}, pixels{std::move(other.pixels)}, domain_lengths{std::move(other.domain_lengths)}, fields{std::move(other.fields)}, F{other.F}, P{other.P}, K{other.K}, // this seems to segfault under clang if it's not a move materials{std::move(other.materials)}, projection{std::move(other.projection)}, is_cell_splitted{std::move(other.is_cell_splitted)} { } /* ---------------------------------------------------------------------- */ template typename CellBase::Material_t & CellBase::add_material(Material_ptr mat) { this->materials.push_back(std::move(mat)); return *this->materials.back(); } /* ---------------------------------------------------------------------- */ template auto CellBase::get_strain_vector() -> Vector_ref { return this->get_strain().eigenvec(); } /* ---------------------------------------------------------------------- */ template auto CellBase::get_stress_vector() const -> ConstVector_ref { return this->get_stress().eigenvec(); } /* ---------------------------------------------------------------------- */ template void CellBase:: set_uniform_strain(const Eigen::Ref & strain) { this->F.get_map() = strain; } /* ---------------------------------------------------------------------- */ template auto CellBase::evaluate_stress() -> ConstVector_ref { if (not this->initialised) { this->initialise(); } for (auto & mat: this->materials) { mat->compute_stresses(this->F, this->P, this->get_formulation()); } return this->P.const_eigenvec(); } /* ---------------------------------------------------------------------- */ template auto CellBase:: evaluate_stress_tangent() -> std::array { if (not this->initialised) { this->initialise(); } constexpr bool create_tangent{true}; this->get_tangent(create_tangent); for (auto & mat: this->materials) { mat->compute_stresses_tangent(this->F, this->P, this->K.value(), this->get_formulation()); } const TangentField_t & k = this->K.value(); return std::array{ this->P.const_eigenvec(), k.const_eigenvec()}; } /* ---------------------------------------------------------------------- */ template auto CellBase:: evaluate_projected_directional_stiffness (Eigen::Ref delF) -> Vector_ref { // the following const_cast should be safe, as long as the // constructed delF_field is const itself const TypedField delF_field ("Proxied raw memory for strain increment", *this->fields, Eigen::Map(const_cast(delF.data()), delF.size()), this->F.get_nb_components()); if (!this->K) { throw std::runtime_error ("currently only implemented for cases where a stiffness matrix " "exists"); } if (delF.size() != this->get_nb_dof()) { std::stringstream err{}; err << "input should be of size ndof = ¶(" << this->subdomain_resolutions << ") × " << DimS << "² = "<< this->get_nb_dof() << " but I got " << delF.size(); throw std::runtime_error(err.str()); } const std::string out_name{"δP; temp output for directional stiffness"}; auto & delP = this->get_managed_field(out_name); auto Kmap{this->K.value().get().get_map()}; auto delPmap{delP.get_map()}; MatrixFieldMap delFmap(delF_field); for (auto && tup: akantu::zip(Kmap, delFmap, delPmap)) { auto & k = std::get<0>(tup); auto & df = std::get<1>(tup); auto & dp = std::get<2>(tup); dp = Matrices::tensmult(k, df); } return Vector_ref(this->project(delP).data(), this->get_nb_dof()); } /* ---------------------------------------------------------------------- */ template std::array CellBase::get_strain_shape() const { return this->projection->get_strain_shape(); } /* ---------------------------------------------------------------------- */ template void CellBase::apply_projection(Eigen::Ref vec) { TypedField field("Proxy for projection", *this->fields, vec, this->F.get_nb_components()); this->projection->apply_projection(field); } /* ---------------------------------------------------------------------- */ template typename CellBase::FullResponse_t CellBase:: evaluate_stress_tangent(StrainField_t & grad) { if (this->initialised == false) { this->initialise(); } //! High level compatibility checks if (grad.size() != this->F.size()) { throw std::runtime_error("Size mismatch"); } constexpr bool create_tangent{true}; this->get_tangent(create_tangent); for (auto & mat: this->materials) { mat->compute_stresses_tangent(grad, this->P, this->K.value(), - this->get_formulation(), this->is_cell_splitted); + this->get_formulation()); } return std::tie(this->P, this->K.value()); } /* ---------------------------------------------------------------------- */ template typename CellBase::StressField_t & CellBase::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 CellBase::Vector_ref CellBase::directional_stiffness_vec(const Eigen::Ref &delF) { if (!this->K) { throw std::runtime_error ("currently only implemented for cases where a stiffness matrix " "exists"); } if (delF.size() != this->get_nb_dof()) { std::stringstream err{}; err << "input should be of size ndof = ¶(" << this->subdomain_resolutions << ") × " << DimS << "² = "<< this->get_nb_dof() << " but I got " << delF.size(); throw std::runtime_error(err.str()); } const std::string out_name{"temp output for directional stiffness"}; const std::string in_name{"temp input for directional stiffness"}; auto & out_tempref = this->get_managed_field(out_name); auto & in_tempref = this->get_managed_field(in_name); Vector_ref(in_tempref.data(), this->get_nb_dof()) = delF; this->directional_stiffness(this->K.value(), in_tempref, out_tempref); return Vector_ref(out_tempref.data(), this->get_nb_dof()); } /* ---------------------------------------------------------------------- */ template Eigen::ArrayXXd CellBase:: directional_stiffness_with_copy (Eigen::Ref delF) { if (!this->K) { throw std::runtime_error ("currently only implemented for cases where a stiffness matrix " "exists"); } const std::string out_name{"temp output for directional stiffness"}; const std::string in_name{"temp input for directional stiffness"}; auto & out_tempref = this->get_managed_field(out_name); auto & in_tempref = this->get_managed_field(in_name); in_tempref.eigen() = delF; this->directional_stiffness(this->K.value(), in_tempref, out_tempref); return out_tempref.eigen(); } /* ---------------------------------------------------------------------- */ template typename CellBase::StressField_t & CellBase::project(StressField_t &field) { this->projection->apply_projection(field); return field; } /* ---------------------------------------------------------------------- */ template typename CellBase::StrainField_t & CellBase::get_strain() { if (this->initialised == false) { this->initialise(); } return this->F; } /* ---------------------------------------------------------------------- */ template const typename CellBase::StressField_t & CellBase::get_stress() const { return this->P; } /* ---------------------------------------------------------------------- */ template const typename CellBase::TangentField_t & CellBase::get_tangent(bool create) { if (!this->K) { if (create) { this->K = make_field("Tangent Stiffness", *this->fields); } else { throw std::runtime_error ("K does not exist"); } } return this->K.value(); } /* ---------------------------------------------------------------------- */ template typename CellBase::StrainField_t & CellBase::get_managed_field(std::string unique_name) { if (!this->fields->check_field_exists(unique_name)) { return make_field(unique_name, *this->fields); } else { return static_cast(this->fields->at(unique_name)); } } /* ---------------------------------------------------------------------- */ template void CellBase::initialise(FFT_PlanFlags flags) { // check that all pixels have been assigned exactly one material this->check_material_coverage(); for (auto && mat: this->materials) { mat->initialise(); } // resize all global fields (strain, stress, etc) this->fields->initialise(this->subdomain_resolutions, this->subdomain_locations); // initialise the projection and compute the fft plan this->projection->initialise(flags); this->initialised = true; } /* ---------------------------------------------------------------------- */ template void CellBase::save_history_variables() { for (auto && mat: this->materials) { mat->save_history_variables(); } } /* ---------------------------------------------------------------------- */ template typename CellBase::iterator CellBase::begin() { return this->pixels.begin(); } /* ---------------------------------------------------------------------- */ template typename CellBase::iterator CellBase::end() { return this->pixels.end(); } /* ---------------------------------------------------------------------- */ template auto CellBase::get_adaptor() -> Adaptor { return Adaptor(*this); } /* ---------------------------------------------------------------------- */ template void CellBase::check_material_coverage() { auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector*> assignments(nb_pixels, nullptr); for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, 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->subdomain_resolutions, this->subdomain_locations, 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 Rcoord_t CellBase::get_pixel_lengths() const{ auto nb_pixels = this->get_domain_resolutions(); auto length_pixels = this->get_domain_lengths(); Rcoord_t ret_val; for (int i = 0 ; i < DimS ; i++ ){ ret_val[i] = length_pixels[i] / nb_pixels[i]; } return ret_val; } template Rcoord_t CellBase::get_pixel_coordinate(Ccoord_t & pixel) const{ auto pixel_length = this->get_pixel_lengths(); Rcoord_t pixel_coordinate; for (int i = 0 ; i < DimS; i++){ pixel_coordinate[i] = pixel_length[i] * (pixel[i] + 0.5); //0.5 is added to shift to the center of the pixel; } // return pixel_length * pixel; return pixel_coordinate; } template bool CellBase::is_inside(Rcoord_t point){ auto length_pixels = this->get_domain_lengths(); Dim_t counter = 1; for(int i = 0; i < DimS ; i++){ if(point[i] < length_pixels[i]) { counter++; } } return counter == DimS; } template bool CellBase::is_inside(Ccoord_t pixel){ auto nb_pixels = this->get_domain_resolutions(); Dim_t counter = 1; for(int i = 0; i < DimS ; i++){ if(pixel[i] < nb_pixels[i]) { counter++; } } return counter == DimS; } template Real CellBase::get_pixel_volume(){ auto pixel_length = get_pixel_lengths(); Real Retval = 1.0; for (int i = 0 ; i < DimS; i++){ Retval*=pixel_length[i]; } return Retval; } template class CellBase; template class CellBase; } // muSpectre diff --git a/src/cell/cell_base.hh b/src/cell/cell_base.hh index 384554b..5f7046e 100644 --- a/src/cell/cell_base.hh +++ b/src/cell/cell_base.hh @@ -1,523 +1,523 @@ /** * @file cell_base.hh * * @author Till Junge * * @date 01 Nov 2017 * * @brief Base class representing a unit cell cell with single * projection operator * * 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 CELL_BASE_H #define CELL_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 "cell/cell_traits.hh" #include #include #include #include #include namespace muSpectre { /** * Cell adaptors implement the matrix-vector multiplication and * allow the system to be used like a sparse matrix in * conjugate-gradient-type solvers */ template class CellAdaptor; /** * Base class for cells that is not templated and therefore can be * in solvers that see cells as runtime-polymorphic objects. This * allows the use of standard * (i.e. spectral-method-implementation-agnostic) solvers, as for * instance the scipy solvers */ class Cell { public: //! sparse matrix emulation using Adaptor = CellAdaptor; //! dynamic vector type for interactions with numpy/scipy/solvers etc. using Vector_t = Eigen::Matrix; //! dynamic matrix type for setting strains using Matrix_t = Eigen::Matrix; //! ref to constant vector using ConstVector_ref = Eigen::Map; //! output vector reference for solvers using Vector_ref = Eigen::Map; //! Default constructor Cell() = default; //! Copy constructor Cell(const Cell &other) = default; //! Move constructor Cell(Cell &&other) = default; //! Destructor virtual ~Cell() = default; //! Copy assignment operator Cell& operator=(const Cell &other) = default; //! Move assignment operator Cell& operator=(Cell &&other) = default; //! for handling double initialisations right bool is_initialised() const {return this->initialised;} //! returns the number of degrees of freedom in the cell virtual Dim_t get_nb_dof() const = 0; //! return the communicator object virtual const Communicator & get_communicator() const = 0; /** * formulation is hard set by the choice of the projection class */ virtual const Formulation & get_formulation() const = 0; /** * returns the material dimension of the problem */ virtual Dim_t get_material_dim() const = 0; /** * returns the number of rows and cols for the strain matrix type * (for full storage, the strain is stored in material_dim × * material_dim matrices, but in symmetriy storage, it is a column * vector) */ virtual std::array get_strain_shape() const = 0; /** * returns a writable map onto the strain field of this cell. This * corresponds to the unknowns in a typical solve cycle. */ virtual Vector_ref get_strain_vector() = 0; /** * returns a read-only map onto the stress field of this * cell. This corresponds to the intermediate (and finally, total) * solution in a typical solve cycle */ virtual ConstVector_ref get_stress_vector() const = 0; /** * evaluates and returns the stress for the currently set strain */ virtual ConstVector_ref evaluate_stress() = 0; /** * evaluates and returns the stress and stiffness for the currently set strain */ virtual std::array evaluate_stress_tangent() = 0; /** * applies the projection operator in-place on the input vector */ virtual void apply_projection(Eigen::Ref vec) = 0; /** * freezes all the history variables of the materials */ virtual void save_history_variables() = 0; /** * evaluates the directional and projected stiffness (this * corresponds to G:K:δF in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that * this operation needs to be implemented with a copy in oder to * be compatible with scipy and EigenCG etc (At the very least, * the copy is only made once) */ virtual Vector_ref evaluate_projected_directional_stiffness (Eigen::Ref delF) = 0; /** * set uniform strain (typically used to initialise problems */ virtual void set_uniform_strain(const Eigen::Ref &) = 0; //! get a sparse matrix view on the cell virtual Adaptor get_adaptor() = 0; protected: bool initialised{false}; //!< to handle double initialisation right private: }; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class CellBase: public Cell { public: using Parent = Cell; using Ccoord = Ccoord_t; //!< cell coordinates type using Rcoord = Rcoord_t; //!< physical coordinates type //! global field collection using FieldCollection_t = GlobalFieldCollection; //! the collection is handled in a `std::unique_ptr` using Collection_ptr = std::unique_ptr; //! polymorphic base material type using Material_t = MaterialBase; //! materials handled through `std::unique_ptr`s using Material_ptr = std::unique_ptr; //! polymorphic base projection type using Projection_t = ProjectionBase; //! projections handled through `std::unique_ptr`s using Projection_ptr = std::unique_ptr; //! expected type for strain fields using StrainField_t = TensorField; //! expected type for stress fields using StressField_t = TensorField; //! expected type for tangent stiffness fields using TangentField_t = TensorField; //! combined stress and tangent field using FullResponse_t = std::tuple; //! iterator type over all cell pixel's using iterator = typename CcoordOps::Pixels::iterator; //! dynamic vector type for interactions with numpy/scipy/solvers etc. using Vector_t = typename Parent::Vector_t; //! ref to constant vector using ConstVector_ref = typename Parent::ConstVector_ref; //! output vector reference for solvers using Vector_ref = typename Parent::Vector_ref; //! sparse matrix emulation using Adaptor = Parent::Adaptor; //! Default constructor CellBase() = delete; //! constructor using sizes and resolution - CellBase(Projection_ptr projection, SplittedCell is_cell_splitted = SplittedCell::no); + CellBase(Projection_ptr projection); //! Copy constructor CellBase(const CellBase &other) = delete; //! Move constructor CellBase(CellBase &&other); //! Destructor virtual ~CellBase() = default; //! Copy assignment operator CellBase& operator=(const CellBase &other) = delete; //! Move assignment operator CellBase& operator=(CellBase &&other) = default; /** * Materials can only be moved. This is to assure exclusive * ownership of any material by this cell */ Material_t & add_material(Material_ptr mat); /** * returns a writable map onto the strain field of this cell. This * corresponds to the unknowns in a typical solve cycle. */ virtual Vector_ref get_strain_vector() override; /** * returns a read-only map onto the stress field of this * cell. This corresponds to the intermediate (and finally, total) * solution in a typical solve cycle */ virtual ConstVector_ref get_stress_vector() const override; /** * evaluates and returns the stress for the currently set strain */ virtual ConstVector_ref evaluate_stress() override; /** * evaluates and returns the stress and stiffness for the currently set strain */ virtual std::array evaluate_stress_tangent() override; /** * evaluates the directional and projected stiffness (this * corresponds to G:K:δF in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that * this operation needs to be implemented with a copy in oder to * be compatible with scipy and EigenCG etc. (At the very least, * the copy is only made once) */ virtual Vector_ref evaluate_projected_directional_stiffness (Eigen::Ref delF) override; //! return the template param DimM (required for polymorphic use of `Cell` Dim_t get_material_dim() const override final {return DimM;} /** * returns the number of rows and cols for the strain matrix type * (for full storage, the strain is stored in material_dim × * material_dim matrices, but in symmetriy storage, it is a column * vector) */ std::array get_strain_shape() const override final; /** * applies the projection operator in-place on the input vector */ void apply_projection(Eigen::Ref vec) override final; /** * set uniform strain (typically used to initialise problems */ void set_uniform_strain(const Eigen::Ref &) override; /** * evaluate all materials */ virtual FullResponse_t evaluate_stress_tangent(StrainField_t & F); /** * evaluate directional stiffness (i.e. G:K:δF or G:K:δε) */ StressField_t & directional_stiffness(const TangentField_t & K, const StrainField_t & delF, StressField_t & delP); /** * vectorized version for eigen solvers, no copy, but only works * when fields have ArrayStore=false */ Vector_ref directional_stiffness_vec(const Eigen::Ref & delF); /** * Evaluate directional stiffness into a temporary array and * return a copy. This is a costly and wasteful interface to * directional_stiffness and should only be used for debugging or * in the python interface */ Eigen::ArrayXXd directional_stiffness_with_copy (Eigen::Ref delF); /** * Convenience function circumventing the neeed to use the * underlying projection */ StressField_t & project(StressField_t & field); //! returns a ref to the cell's strain field StrainField_t & get_strain(); //! returns a ref to the cell's stress field const StressField_t & get_stress() const; //! returns a ref to the cell's tangent stiffness field const TangentField_t & get_tangent(bool create = false); //! returns a ref to a temporary field managed by the cell StrainField_t & get_managed_field(std::string unique_name); /** * 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); /** * for materials with state variables, these typically need to be * saved/updated an the end of each load increment, this function * calls this update for each material in the cell */ void save_history_variables() override final; iterator begin(); //!< iterator to the first pixel iterator end(); //!< iterator past the last pixel //! number of pixels in the cell size_t size() const {return pixels.size();} //! return the subdomain resolutions of the cell const Ccoord & get_subdomain_resolutions() const { return this->subdomain_resolutions;} //! return the subdomain locations of the cell const Ccoord & get_subdomain_locations() const { return this->subdomain_locations;} //! return the domain resolutions of the cell const Ccoord & get_domain_resolutions() const { return this->domain_resolutions;} //! return the sizes of the cell const Rcoord & get_domain_lengths() const {return this->domain_lengths;} //! return the real space coordinates of the center of pixel Rcoord get_pixel_coordinate(Ccoord & pixel) const; //! return the physical dimension of a pixel Rcoord get_pixel_lengths() const; //! check if the pixel is inside of the cell bool is_inside(Rcoord point); //! check if the point is inside of the cell bool is_inside(Ccoord pixel); //calculate the volume of each pixel: Real get_pixel_volume(); /** * formulation is hard set by the choice of the projection class */ const Formulation & get_formulation() const override final { return this->projection->get_formulation();} /** * get a reference to the projection object. should only be * required for debugging */ Eigen::Map get_projection() { return this->projection->get_operator();} //! returns the spatial size constexpr static Dim_t get_sdim() {return DimS;}; //! return a sparse matrix adaptor to the cell Adaptor get_adaptor() override; //! returns the number of degrees of freedom in the cell Dim_t get_nb_dof() const override {return this->size()*ipow(DimS, 2);}; //! return the communicator object virtual const Communicator & get_communicator() const override { return this->projection->get_communicator(); } protected: //! make sure that every pixel is assigned to one and only one material virtual void check_material_coverage(); const Ccoord & subdomain_resolutions; //!< the cell's subdomain resolutions const Ccoord & subdomain_locations; //!< the cell's subdomain resolutions const Ccoord & domain_resolutions; //!< the cell's domain resolutions CcoordOps::Pixels pixels; //!< helper to iterate over the pixels const Rcoord & domain_lengths; //!< the cell's lengths Collection_ptr fields; //!< handle for the global fields of the cell StrainField_t & F; //!< ref to strain field StressField_t & P; //!< ref to stress field //! Tangent field might not even be required; so this is an //! optional ref_wrapper instead of a ref optional> K{}; //! container of the materials present in the cell std::vector materials{}; Projection_ptr projection; //!< handle for the projection operator bool initialised{false}; //!< to handle double initialisation right - SplittedCell is_cell_splitted ; + SplittedCell is_cell_splitted {SplittedCell::no}; private: }; /** * lightweight resource handle wrapping a `muSpectre::Cell` or * a subclass thereof into `Eigen::EigenBase`, so it can be * interpreted as a sparse matrix by Eigen solvers */ template class CellAdaptor: public Eigen::EigenBase> { public: using Scalar = double; //!< sparse matrix traits using RealScalar = double; //!< sparse matrix traits using StorageIndex = int; //!< sparse matrix traits enum { ColsAtCompileTime = Eigen::Dynamic, MaxColsAtCompileTime = Eigen::Dynamic, RowsAtCompileTime = Eigen::Dynamic, MaxRowsAtCompileTime = Eigen::Dynamic, IsRowMajor = false }; //! constructor CellAdaptor(Cell & cell):cell{cell}{} //!returns the number of logical rows Eigen::Index rows() const {return this->cell.get_nb_dof();} //!returns the number of logical columns Eigen::Index cols() const {return this->rows();} //! implementation of the evaluation template Eigen::Product operator*(const Eigen::MatrixBase& x) const { return Eigen::Product (*this, x.derived()); } Cell & cell; //!< ref to the cell }; } // muSpectre namespace Eigen { namespace internal { //! Implementation of `muSpectre::CellAdaptor` * `Eigen::DenseVector` through a //! specialization of `Eigen::internal::generic_product_impl`: template struct generic_product_impl // GEMV stands for matrix-vector : generic_product_impl_base > { //! undocumented typedef typename Product::Scalar Scalar; //! undocumented template static void scaleAndAddTo(Dest& dst, const CellAdaptor& lhs, const Rhs& rhs, const Scalar& /*alpha*/) { // This method should implement "dst += alpha * lhs * rhs" inplace, // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it. // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs, dst.noalias() += const_cast(lhs).cell. evaluate_projected_directional_stiffness(rhs); } }; } } #endif /* CELL_BASE_H */ diff --git a/src/cell/cell_split.cc b/src/cell/cell_split.cc index b25254d..8deaee3 100644 --- a/src/cell/cell_split.cc +++ b/src/cell/cell_split.cc @@ -1,282 +1,315 @@ /** * @file cell_split.cc * * @author Ali Falsafi * * @date 19 Apr 2018 * * @brief Implementation for cell base class * * 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 "common/ccoord_operations.hh" #include "common/field.hh" #include "common/utilities.hh" #include "materials/material_base.hh" #include "fft/projection_base.hh" #include "cell/cell_traits.hh" #include "cell/cell_base.hh" #include "cell/cell_split.hh" #include #include #include #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template CellSplit::CellSplit(Projection_ptr projection, SplittedCell is_cell_splitted ) - :Parent(std::move(projection), is_cell_splitted) - {} + :Parent(std::move(projection)) + { + this->is_cell_splitted = is_cell_splitted; + } /* ---------------------------------------------------------------------- */ template std::vector CellSplit::get_assigned_ratios(){ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector pixel_assigned_ratios(nb_pixels, 0.0); for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel); } } return pixel_assigned_ratios; } /* ---------------------------------------------------------------------- */ template std::vector CellSplit::get_unassigned_ratios_incomplete_pixels(){ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector pixel_assigned_ratios(nb_pixels, 0.0); std::vector pixel_assigned_ratios_incomplete_pixels; for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel); } } for (auto && ratio:pixel_assigned_ratios){ if (ratio<1){ pixel_assigned_ratios_incomplete_pixels.push_back(1.0-ratio); } } return pixel_assigned_ratios_incomplete_pixels; } /* ---------------------------------------------------------------------- */ template auto CellSplit::make_incomplete_pixels() -> IncompletePixels{ return IncompletePixels(*this); } /* ---------------------------------------------------------------------- */ template std::vector CellSplit::get_index_incomplete_pixels(){ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector pixel_assigned_ratios(nb_pixels, 0.0); std::vector index_unassigned_pixels; for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel); } } for (auto && tup: akantu::enumerate(this->pixels)){ auto && i{std::get<0> (tup)}; if (pixel_assigned_ratios[i] < 1){ index_unassigned_pixels.push_back(i); } } return index_unassigned_pixels; } /* ---------------------------------------------------------------------- */ template std::vector> CellSplit::get_unassigned_pixels(){ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector pixel_assigned_ratios(nb_pixels, 0.0); std::vector> unassigned_pixels; for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel); } } for (auto && tup: akantu::enumerate(this->pixels)){ auto && index{std::get<0> (tup)}; auto && pix{std::get<1> (tup)}; if (pixel_assigned_ratios[index] < 1){ unassigned_pixels.push_back(pix); } } return unassigned_pixels; } /* ---------------------------------------------------------------------- */ /*template template void CellSplit::complete_material_assignment(MaterialType& material){ std::vector pixel_assigned_ratio(this->get_assigned_ratios()); for (auto && tup: akantu::enumerate(*this)) { auto && pixel = std::get<1>(tup); auto iterator = std::get<0>(tup); if (pixel_assigned_ratio[iterator] < 1.0){ material.add_pixel_split(pixel, 1.0-pixel_assigned_ratio[iterator]); } } }*/ template void CellSplit::check_material_coverage(){ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector pixel_assigned_ratio(nb_pixels, 0.0); for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); pixel_assigned_ratio.at(index) += mat->get_assigned_ratio(pixel); } } std::vector> over_assigned_pixels; std::vector> under_assigned_pixels; for (size_t i = 0; i < nb_pixels; ++i) { if (pixel_assigned_ratio.at(i) > 1.0) { over_assigned_pixels.push_back(CcoordOps::get_ccoord(this->subdomain_resolutions, this->subdomain_locations, i)); }else if (pixel_assigned_ratio.at(i) < 1.0) { under_assigned_pixels.push_back(CcoordOps::get_ccoord(this->subdomain_resolutions, this->subdomain_locations, i)); } } if (over_assigned_pixels.size() != 0) { std::stringstream err {}; err << "Execesive material is assigned to the following pixels: "; for (auto & pixel: over_assigned_pixels) { err << pixel << ", "; } err << "and that cannot be handled"; throw std::runtime_error(err.str()); } if (under_assigned_pixels.size() != 0) { std::stringstream err {}; err << "Insufficient material is assigned to the following pixels: "; for (auto & pixel: under_assigned_pixels) { err << pixel << ", "; } err << "and that cannot be handled"; throw std::runtime_error(err.str()); } } /* ---------------------------------------------------------------------- */ //this piece of code handles the evaluation of stress an dtangent matrix //in case the cells have materials in which pixels are partially composed of //diffferent materials. template typename CellSplit::FullResponse_t CellSplit::evaluate_stress_tangent(StrainField_t & F){ - return evaluate_split_stress_tangent(F); - } - /* ---------------------------------------------------------------------- */ - template - typename CellSplit::FullResponse_t - CellSplit::evaluate_split_stress_tangent(StrainField_t & F){ if (this->initialised == false) { this->initialise(); } //! High level compatibility checks if (F.size() != this->F.size()) { throw std::runtime_error("Size mismatch"); } constexpr bool create_tangent{true}; this->get_tangent(create_tangent); // Here we should first set P and K matrixes equal to zeros first because we want to add up contribution //of the partial influence of different materials assigend to each pixel. Therefore, this values should be // initiialised as zero filled tensors this->set_p_k_zero(); for (auto & mat: this->materials) { mat->compute_stresses_tangent(F, this->P, this->K.value(), this->get_formulation(), this->is_cell_splitted); } return std::tie(this->P, this->K.value()); + } + /* ---------------------------------------------------------------------- */ + template + auto CellSplit:: + evaluate_stress() -> ConstVector_ref { + if (not this->initialised) { + this->initialise(); + } + this->P.set_zero(); + for (auto & mat: this->materials) { + mat->compute_stresses(this->F, this->P, this->get_formulation()); + } + + return this->P.const_eigenvec(); + } +/* ---------------------------------------------------------------------- */ + template + auto CellSplit:: + evaluate_stress_tangent() -> std::array { + if (not this->initialised) { + this->initialise(); + } + + constexpr bool create_tangent{true}; + this->get_tangent(create_tangent); + + // Here we should first set P and K matrixes equal to zeros first because we want to add up contribution + //of the partial influence of different materials assigend to each pixel. Therefore, this values should be + // initiialised as zero filled tensors + this->set_p_k_zero(); + for (auto & mat: this->materials) { + mat->compute_stresses_tangent(this->F, this->P, this->K.value(), + this->get_formulation(), this->is_cell_splitted); + } + const TangentField_t & k = this->K.value(); + return std::array{ + this->P.const_eigenvec(), k.const_eigenvec()}; + } /* ---------------------------------------------------------------------- */ template void CellSplit::set_p_k_zero(){ this->P.set_zero(); this->K.value().get().set_zero(); } - /* ---------------------------------------------------------------------- */ template CellSplit::IncompletePixels::IncompletePixels(CellSplit& cell) :cell(cell),incomplete_assigned_ratios(cell.get_unassigned_ratios_incomplete_pixels()), index_incomplete_pixels(cell.get_index_incomplete_pixels()) {} /* ---------------------------------------------------------------------- */ template CellSplit::IncompletePixels::iterator::iterator (const IncompletePixels & pixels, bool begin) : incomplete_pixels(pixels), index{begin? 0:this->incomplete_pixels.index_incomplete_pixels.size()} {} /* ---------------------------------------------------------------------- */ template auto CellSplit::IncompletePixels::iterator::operator++()->iterator&{ this->index++; return *this; } /* ---------------------------------------------------------------------- */ template auto CellSplit::IncompletePixels::iterator::operator!=(const iterator & other)->bool{ // return (this->incomplete_pixels.index_incomplete_pixels[index] != // other.incomplete_pixels.index_incomplete_pixels[index]); return (this->index != other.index); } /* ---------------------------------------------------------------------- */ template auto CellSplit::IncompletePixels::iterator::operator*()->value_type const{ auto ccoord = CcoordOps::get_ccoord(this->incomplete_pixels.cell.get_domain_resolutions(), this->incomplete_pixels.cell.get_subdomain_locations(), this->incomplete_pixels.index_incomplete_pixels[index]); auto ratio = this->incomplete_pixels.incomplete_assigned_ratios[index]; return std::make_tuple(ccoord, ratio); } /* ---------------------------------------------------------------------- */ template class CellSplit; template class CellSplit; /* ---------------------------------------------------------------------- */ } //muspectre diff --git a/src/cell/cell_split.hh b/src/cell/cell_split.hh index 20820b6..808ac8d 100644 --- a/src/cell/cell_split.hh +++ b/src/cell/cell_split.hh @@ -1,193 +1,197 @@ /** * @file cell_split.hh * * @author Ali Falsafi * * @date 19 Apr 2018 * * @brief Base class representing a unit cell able to handle * split material assignments * * 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 CELL_SPLIT_H #define CELL_SPLIT_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 "cell/cell_traits.hh" #include "cell/cell_base.hh" #include "common/intersection_octree.hh" #include #include #include #include namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) /* This class handles the cells that has splitly assigned material to their pixels */ template class CellSplit: public CellBase { friend class CellBase; public: using Parent = CellBase; //!< base class //! global field collection using FieldCollection_t = GlobalFieldCollection; using Projection_t = ProjectionBase; //! projections handled through `std::unique_ptr`s using Projection_ptr = std::unique_ptr; using StrainField_t = TensorField; //! expected type for stress fields using StressField_t = TensorField; //! expected type for tangent stiffness fields using TangentField_t = TensorField; //! combined stress and tangent field using FullResponse_t = std::tuple; + //! ref to constant vector + using ConstVector_ref = typename Parent::ConstVector_ref; + //! Default constructor CellSplit() = delete; //! constructor using sizes and resolution CellSplit(Projection_ptr projection, SplittedCell is_cell_splitted = SplittedCell::yes); //! Copy constructor CellSplit(const CellSplit &other) = delete; //! Move constructor CellSplit(CellSplit &&other) = default; //! Destructor virtual ~CellSplit() = default; //! Copy assignment operator CellSplit& operator=(const CellSplit &other) = delete; //! Move assignment operator CellSplit& operator=(CellSplit &&other) = default; //completes the assignmnet of material with a specific material so all the under-assigned pixels would be assigned to a material void complete_material_assignment(MaterialBase& material); //returns the assigend ratios to each pixel std::vector get_assigned_ratios(); // //template void make_automatic_precipitate(std::vector> preciptiate_vertices, MaterialBase& material); // std::vector get_unassigned_ratios_incomplete_pixels(); std::vector get_index_incomplete_pixels(); std::vector> get_unassigned_pixels(); class IncompletePixels{ public: //! constructor IncompletePixels(CellSplit& cell); //! copy constructor IncompletePixels(const IncompletePixels & other) = default; //! move constructor IncompletePixels(IncompletePixels & other) = default; //Deconstructor virtual ~IncompletePixels() = default; //! iterator type over all incompletetedly assigned pixel's class iterator{ public: using value_type = std::tuple, Real>; //!< stl conformance using const_value_type = const value_type; //!< stl conformance using pointer = value_type*; //!< stl conformance using difference_type = std::ptrdiff_t; //!< stl conformance using iterator_category = std::forward_iterator_tag;//!value_type const; //! pre-increment auto operator++()->iterator&; //! inequality auto operator!=(const iterator & other)->bool; //! equality inline bool operator==(const iterator & other) const; private: const IncompletePixels& incomplete_pixels; size_t index; }; inline iterator begin() const {return iterator(*this);} //! stl conformance inline iterator end() const {return iterator(*this, false);} //! stl conformance inline size_t size() const {return CcoordOps::get_size(this->cell.get_domain_resolutions());} private: CellSplit& cell; std::vector incomplete_assigned_ratios; std::vector index_incomplete_pixels; }; auto make_incomplete_pixels() -> IncompletePixels; protected: void check_material_coverage() override final; void set_p_k_zero(); //full resppnse is consisted of the stresses and tangent matrix FullResponse_t evaluate_stress_tangent(StrainField_t & F) override final; - FullResponse_t evaluate_split_stress_tangent(StrainField_t & F); + std::array evaluate_stress_tangent() override final; + ConstVector_ref evaluate_stress() override final; private: }; /* ---------------------------------------------------------------------- */ template void CellSplit::make_automatic_precipitate(std::vector> precipitate_vertices, MaterialBase& material){ RootNode precipitate (*this, precipitate_vertices); auto precipitate_intersects= precipitate.get_intersected_pixels(); auto precipitate_intersection_ratios = precipitate.get_intersection_ratios(); for (auto tup:akantu::zip(precipitate_intersects, precipitate_intersection_ratios)){ auto pix = std::get<0> (tup); auto ratio = std::get<1> (tup); material.add_pixel_split(pix,ratio); } } /* ---------------------------------------------------------------------- */ template void CellSplit::complete_material_assignment(MaterialBase& material){ std::vector pixel_assigned_ratio(this->get_assigned_ratios()); for (auto && tup: akantu::enumerate(*this)) { auto && pixel = std::get<1>(tup); auto iterator = std::get<0>(tup); if (pixel_assigned_ratio[iterator] < 1.0){ material.add_pixel_split(pixel, 1.0-pixel_assigned_ratio[iterator]); } } } } //muspectre #endif /* CELL_SPLIT_H */ diff --git a/src/materials/CMakeLists.txt b/src/materials/CMakeLists.txt index 722e08c..506bd0f 100644 --- a/src/materials/CMakeLists.txt +++ b/src/materials/CMakeLists.txt @@ -1,41 +1,43 @@ # ============================================================================= # file CMakeLists.txt # # @author Till Junge # # @date 08 Jan 2018 # # @brief configuration for material laws # # @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 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. # ============================================================================= set (materials_SRC ${CMAKE_CURRENT_SOURCE_DIR}/material_base.cc + ${CMAKE_CURRENT_SOURCE_DIR}/material_anisotropic.cc + ${CMAKE_CURRENT_SOURCE_DIR}/material_orthotropic.cc ${CMAKE_CURRENT_SOURCE_DIR}/material_linear_elastic1.cc ${CMAKE_CURRENT_SOURCE_DIR}/material_linear_elastic2.cc ${CMAKE_CURRENT_SOURCE_DIR}/material_linear_elastic3.cc ${CMAKE_CURRENT_SOURCE_DIR}/material_linear_elastic4.cc ${CMAKE_CURRENT_SOURCE_DIR}/material_hyper_elasto_plastic1.cc ) target_sources(muSpectre PRIVATE ${materials_SRC}) diff --git a/src/materials/material_anisotropic.cc b/src/materials/material_anisotropic.cc new file mode 100644 index 0000000..4595d31 --- /dev/null +++ b/src/materials/material_anisotropic.cc @@ -0,0 +1,104 @@ +/** + * @file material_anisotropic.cc + * + * @author Ali Falsafi + * + * @date 09 Jul 2018 + * + * @brief Base class for materials (constitutive models) + * + * 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 "material_base.hh" +#include "common/common.hh" +#include "material_anisotropic.hh" +#include "common/T4_map_proxy.hh" +#include "common/tensor_algebra.hh" +#include "common/eigen_tools.hh" + +namespace muSpectre { + /* ---------------------------------------------------------------------- */ + template + MaterialAnisotropic::MaterialAnisotropic(std::string name, + std::vector input_c) + :Parent(name), C{c_maker(input_c)} + {}; + + + + /* ---------------------------------------------------------------------- */ + template + auto MaterialAnisotropic::c_maker(std::vector input) + -> Stiffness_t + { + //Stiffness_t Ctmp; + //T4MatMap C4{Ctmp.data()}; + Stiffness_t C4; + if (DimS == twoD){ + try{ + if (input.size() != 6){ + throw 1; + } + } + catch(int x) { + std::cout<<"Number of the entries should be 6"<; +template class MaterialAnisotropic; +template class MaterialAnisotropic; + +}//muSpectre diff --git a/src/materials/material_anisotropic.hh b/src/materials/material_anisotropic.hh new file mode 100644 index 0000000..709f3f9 --- /dev/null +++ b/src/materials/material_anisotropic.hh @@ -0,0 +1,159 @@ +/** + * @file material_anisotropic.hh + * + * @author Ali Falsafi + * + * @date 9 Jul 2018 + * + * @brief Base class for materials (constitutive models) + * + * 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_ANISOTROPIC_H +#define MATERIAL_ANISOTROPIC_H + +#include "materials/material_base.hh" +#include "materials/material_muSpectre_base.hh" +#include "common/common.hh" +#include "common/T4_map_proxy.hh" + +namespace muSpectre { + template + class MaterialAnisotropic; + + //traits for anisotropic material + template + struct MaterialMuSpectre_traits> { + using Parent = MaterialMuSpectre_traits;//!< base for elasticity + + //! global field collection + using GFieldCollection_t = typename + MaterialBase::GFieldCollection_t; + + //! expected map type for strain fields + using StrainMap_t = MatrixFieldMap; + //! expected map type for stress fields + using StressMap_t = MatrixFieldMap; + //! expected map type for tangent stiffness fields + using TangentMap_t = T4MatrixFieldMap; + + //! declare what type of strain measure your law takes as input + constexpr static auto strain_measure{StrainMeasure::GreenLagrange}; + //! declare what type of stress measure your law yields as output + constexpr static auto stress_measure{StressMeasure::PK2}; + + //! anisotropicity without internal variables + using InternalVariables = std::tuple<>; + }; + + /** + * Material implementation for anisotropic constitutive law + */ + template + class MaterialAnisotropic: + public MaterialMuSpectre< MaterialAnisotropic, DimS, DimM> + { + + //! base class + using Parent = MaterialMuSpectre; + + using Stiffness_t =T4Mat; + + //! traits of this material + using traits = MaterialMuSpectre_traits; + + //! this law does not have any internal variables + using InternalVariables = typename traits::InternalVariables; + + public: + using parent = MaterialMuSpectre + , DimS, DimM>; + + //! global field collection + using GFieldCollection_t = typename + MaterialBase::GFieldCollection_t; + + //! expected map type for tangent stiffness fields + using Tangent_t = T4MatrixFieldMap; + + //! Default constructor + MaterialAnisotropic() = delete; + + MaterialAnisotropic(std::string name, std::vector input_c); + + //! Copy constructor + MaterialAnisotropic(const MaterialAnisotropic &other) = delete; + + //! Move constructor + MaterialAnisotropic(MaterialAnisotropic &&other) = delete; + + //! Destructor + virtual ~MaterialAnisotropic() = default; + + template + inline decltype(auto) evaluate_stress(s_t && E); + + /** + * evaluates both second Piola-Kirchhoff stress and stiffness given + * the Green-Lagrange strain (or Cauchy stress and stiffness if + * called with a small strain tensor) and the local stiffness tensor. + */ + template + inline decltype(auto) + evaluate_stress_tangent(s_t && E); + + /** + * return the empty internals tuple + */ + InternalVariables & get_internals() { + return this->internal_variables;}; + + // takes the elements of the C and makes it: + virtual auto c_maker(std::vector input) -> Stiffness_t; + + protected: + + Stiffness_t C; //!< stiffness tensor + + //! empty tuple + InternalVariables internal_variables{}; + }; + + + /* ---------------------------------------------------------------------- */ + template + template + decltype(auto) + MaterialAnisotropic:: + evaluate_stress(s_t && E) { + return Matrices::tensmult(this->C, E); + } + + /* ---------------------------------------------------------------------- */ + template + template + decltype(auto) + MaterialAnisotropic:: + evaluate_stress_tangent(s_t && E) { + return std::make_tuple + (evaluate_stress(E), this->C); + } + +}//muSpectre + +#endif diff --git a/src/materials/material_base.cc b/src/materials/material_base.cc index f26f7de..f96bf03 100644 --- a/src/materials/material_base.cc +++ b/src/materials/material_base.cc @@ -1,92 +1,94 @@ /** * @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 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 "materials/material_base.hh" namespace muSpectre { //----------------------------------------------------------------------------// template MaterialBase::MaterialBase(std::string name) :name(name),assigned_ratio{make_field("assigned ratio", this->internal_fields)}, assigned_ratio_mapped{assigned_ratio}{ 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 < Dim_t DimS, Dim_t DimM> void MaterialBase< DimS, DimM>:: add_pixel_split(const Ccoord &local_ccoord, Real ratio){ - this->internal_fields.add_pixel_split(local_ccoord, ratio); + auto & this_mat = static_cast(*this); + this_mat.add_pixel(local_ccoord); + this->assigned_ratio.push_back(ratio); } /* ---------------------------------------------------------------------- */ template void MaterialBase::compute_stresses(const Field_t & F, Field_t & P, Formulation form, SplittedCell is_cell_splitted) { this->compute_stresses(StrainField_t::check_ref(F), StressField_t::check_ref(P), form, is_cell_splitted); } template Real MaterialBase::get_assigned_ratio(Ccoord pixel){ return this->assigned_ratio_mapped[pixel]; } /* ---------------------------------------------------------------------- */ template void MaterialBase::compute_stresses_tangent(const Field_t & F, Field_t & P, Field_t & K, Formulation form, SplittedCell is_cell_splitted) { this->compute_stresses_tangent(StrainField_t::check_ref(F), StressField_t::check_ref(P), TangentField_t::check_ref(K), form, is_cell_splitted); } template class MaterialBase<2, 2>; template class MaterialBase<2, 3>; template class MaterialBase<3, 3>; } // muSpectre diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh index 86bfb5f..abf79e7 100644 --- a/src/materials/material_base.hh +++ b/src/materials/material_base.hh @@ -1,167 +1,167 @@ /** * @file material_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials (constitutive models) * * 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_BASE_H #define MATERIAL_BASE_H #include "common/common.hh" #include "common/field.hh" #include "common/field_collection.hh" #include namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) and /** * @a DimM is the material dimension (i.e., the dimension of constitutive * law; even for e.g. two-dimensional problems the constitutive law could * live in three-dimensional space for e.g. plane strain or stress problems) */ template class MaterialBase { public: //! typedefs for data handled by this interface //! global field collection for cell-wide fields, like stress, strain, etc using GFieldCollection_t = GlobalFieldCollection; //! field collection for internal variables, such as eigen-strains, //! plastic strains, damage variables, etc, but also for managing which //! pixels the material is responsible for using MFieldCollection_t = LocalFieldCollection; using iterator = typename MFieldCollection_t::iterator; //!< pixel iterator //! polymorphic base class for fields only to be used for debugging using Field_t = internal::FieldBase; //! Full type for stress fields using StressField_t = TensorField; //! Full type for strain fields using StrainField_t = StressField_t; //! Full type for tangent stiffness fields fields using TangentField_t = TensorField; using Ccoord = Ccoord_t; //!< cell coordinates type //! Default constructor MaterialBase() = delete; //! Construct by name MaterialBase(std::string name); //! Copy constructor MaterialBase(const MaterialBase &other) = delete; //! Move constructor MaterialBase(MaterialBase &&other) = delete; //! Destructor virtual ~MaterialBase() = default; //! Copy assignment operator MaterialBase& operator=(const MaterialBase &other) = delete; //! Move assignment operator MaterialBase& operator=(MaterialBase &&other) = delete; /** * take responsibility for a pixel identified by its cell coordinates * WARNING: this won't work for materials with additional info per pixel * (as, e.g. for eigenstrain), we need to pass more parameters. Materials - * of this tye need to overload add_pixel + * of this type need to overload add_pixel */ virtual void add_pixel(const Ccoord & ccooord); virtual void add_pixel_split (const Ccoord &local_ccoord, Real ratio); //! allocate memory, etc, but also: wipe history variables! virtual void initialise() = 0; /** * for materials with state variables, these typically need to be * saved/updated an the end of each load increment, the virtual * base implementation does nothing, but materials with history * variables need to implement this */ virtual void save_history_variables() {}; //! return the material's name const std::string & get_name() const; //! spatial dimension for static inheritance constexpr static Dim_t sdim() {return DimS;} //! material dimension for static inheritance constexpr static Dim_t mdim() {return DimM;} //! computes stress virtual void compute_stresses(const StrainField_t & F, StressField_t & P, Formulation form, SplittedCell is_cell_splitted) = 0; /** * Convenience function to compute stresses, mostly for debugging and * testing. Has runtime-cost associated with compatibility-checking and * conversion of the Field_t arguments that can be avoided by using the * version with strongly typed field references */ void compute_stresses(const Field_t & F, Field_t & P, Formulation form, SplittedCell is_cell_splitted = SplittedCell::no); //! computes stress and tangent moduli virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form, SplittedCell is_cell_splitted) = 0; /** * Convenience function to compute stresses and tangent moduli, mostly for * debugging and testing. Has runtime-cost associated with * compatibility-checking and conversion of the Field_t arguments that can * be avoided by using the version with strongly typed field references */ void compute_stresses_tangent(const Field_t & F, Field_t & P, Field_t & K, Formulation form, SplittedCell is_cell_splitted = SplittedCell::no); // this function return the ratio of which the input pixel is consisted of this material Real get_assigned_ratio(Ccoord pixel); //! iterator to first pixel handled by this material inline iterator begin() {return this->internal_fields.begin();} //! iterator past the last pixel handled by this material inline iterator end() {return this->internal_fields.end();} //! number of pixels assigned to this material inline size_t size() const {return this->internal_fields.size();} protected: const std::string name; //!< material's name (for output and debugging) MFieldCollection_t internal_fields{};//!< storage for internal variables using ScalarField_t = ScalarField, Real> ; ScalarField_t & assigned_ratio; //!< field holding the assigning ratio of the material using ScalarFieldMap_t = ScalarFieldMap, Real, true>; ScalarFieldMap_t assigned_ratio_mapped; private: }; } // muSpectre #endif /* MATERIAL_BASE_H */ diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index 2eb3036..8a34ba8 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,911 +1,912 @@ /** * @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. * * 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 CellBase; /** * material traits are used by `muSpectre::MaterialMuSpectre` to * break the circular dependence created by the curiously recurring * template parameter. These traits must define * - these `muSpectre::FieldMap`s: * - `StrainMap_t`: typically a `muSpectre::MatrixFieldMap` for a * constant second-order `muSpectre::TensorField` * - `StressMap_t`: typically a `muSpectre::MatrixFieldMap` for a * writable secord-order `muSpectre::TensorField` * - `TangentMap_t`: typically a `muSpectre::T4MatrixFieldMap` for a * writable fourth-order `muSpectre::TensorField` * - `strain_measure`: the expected strain type (will be replaced by the * small-strain tensor ε * `muspectre::StrainMeasure::Infinitesimal` in small * strain computations) * - `stress_measure`: the measure of the returned stress. Is used by * `muspectre::MaterialMuSpectre` to transform it into * Cauchy stress (`muspectre::StressMeasure::Cauchy`) in * small-strain computations and into first * Piola-Kirchhoff stress `muspectre::StressMeasure::PK1` * in finite-strain computations * - `InternalVariables`: a tuple of `muSpectre::FieldMap`s containing * internal variables */ template struct MaterialMuSpectre_traits { }; template class MaterialMuSpectre; /** * Base class for most convenient implementation of materials */ template class MaterialMuSpectre: public MaterialBase { public: /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = MatTB::NeedTangent; using Parent = MaterialBase; //!< base class //! global field collection using GFieldCollection_t = typename Parent::GFieldCollection_t; //! expected type for stress fields using StressField_t = typename Parent::StressField_t; //! expected type for strain fields using StrainField_t = typename Parent::StrainField_t; //! expected type for tangent stiffness fields using TangentField_t = typename Parent::TangentField_t; using Ccoord = Ccoord_t; //!< cell coordinates type //! traits for the CRTP subclass 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) = delete; //! Destructor virtual ~MaterialMuSpectre() = default; //! Factory template static Material & make(CellBase & cell, ConstructorArgs &&... args); //! Copy assignment operator MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete; //! Move assignment operator MaterialMuSpectre& operator=(MaterialMuSpectre &&other) = delete; //! allocate memory, etc virtual void initialise() override; // this fucntion is speicalized to assign partial material to a pixel //void add_pixel_split(const Ccoord & local_ccoord, Real ratio = 1.0); template - inline void add_pixel_split(const Ccoord_t & pixel, + void add_pixel_split(const Ccoord_t & pixel, Real ratio, InternalArgs ... args); - inline void add_pixel_split(const Ccoord_t & pixel, - Real ratio) override; + void add_pixel_split(const Ccoord_t & pixel, + Real ratio = 1.0) override; - - void add_split_pixels_precipitate(std::vector> intersected_precipitate, + //add pixels intersecting to material to the material + void add_split_pixels_precipitate(std::vector> intersected_pixles, std::vector intersection_ratios); + //! computes stress using Parent::compute_stresses; using Parent::compute_stresses_tangent; virtual void compute_stresses(const StrainField_t & F, StressField_t & P, Formulation form, SplittedCell is_cell_splitted) override; //! computes stress and tangent modulus virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form, SplittedCell is_cell_splitted)override; 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 + * inheriting classes need to overload this function */ typename traits::InternalVariables& get_internals() { return static_cast(*this).get_internals();} bool is_initialised{false}; //!< to handle double initialisation right private: }; /* ---------------------------------------------------------------------- */ template MaterialMuSpectre:: MaterialMuSpectre(std::string name) :Parent(name){ using stress_compatible = typename traits::StressMap_t:: template is_compatible; using strain_compatible = typename traits::StrainMap_t:: template is_compatible; using tangent_compatible = typename traits::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(CellBase & cell, ConstructorArgs && ... args) { auto mat = std::make_unique(args...); auto & mat_ref = *mat; cell.add_material(std::move(mat)); return mat_ref; } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: add_pixel_split(const Ccoord &local_ccoord, Real ratio) { auto & this_mat = static_cast(*this); this_mat.add_pixel(local_ccoord); this->assigned_ratio.push_back(ratio); } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre:: add_pixel_split(const Ccoord_t & pixel, Real ratio, InternalArgs ... args){ auto & this_mat = static_cast(*this); this_mat.add_pixel(pixel, args...); this->assigned_ratio.push_back(ratio); } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: - add_split_pixels_precipitate(std::vector> intersected_precipitate, + add_split_pixels_precipitate(std::vector> intersected_pixels, std::vector intersection_ratios){ //assign precipitate materials: - for (auto && tup: akantu::zip(intersected_precipitate, intersection_ratios) ){ + for (auto && tup: akantu::zip(intersected_pixels, intersection_ratios) ){ auto pix = std::get<0> (tup); auto ratio = std::get<1> (tup); this->add_pixel_split(pix,ratio); } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: initialise() { 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, SplittedCell is_cell_splitted) { switch (form) { case Formulation::finite_strain: { switch (is_cell_splitted){ case (SplittedCell::no): { this->compute_stresses_worker(F, P); break; } case (SplittedCell::yes): { this->compute_stresses_worker(F, P); break; } } break; } case Formulation::small_strain: { switch (is_cell_splitted){ case (SplittedCell::no): { this->compute_stresses_worker(F, P); break; } case (SplittedCell::yes): { this->compute_stresses_worker(F, P); break; } } 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, SplittedCell is_cell_splitted) { switch (form) { case Formulation::finite_strain: { switch (is_cell_splitted){ case (SplittedCell::no): { this->compute_stresses_worker(F, P, K); break; } case (SplittedCell::yes): { this->compute_stresses_worker(F, P, K); break; } } break; } case Formulation::small_strain: { switch (is_cell_splitted){ case (SplittedCell::no): { this->compute_stresses_worker(F, P, K); break; } case (SplittedCell::yes): { this->compute_stresses_worker(F, P, K); break; } } 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, const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); // return value contains a tuple of rvalue_refs to both stress and tangent moduli // for the case that cells do not contain any splitt cells auto non_split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables](){ Stresses = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); }; // for the case that cells contain splitt cells auto split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables](){ auto stress_tgt_contributions = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); auto && stress = std::get<0>(Stresses); auto && tangent = std::get<1>(Stresses); stress += ratio * std::get<0>(stress_tgt_contributions); tangent += ratio * std::get<1>(stress_tgt_contributions); }; switch (is_cell_splitted){ case SplittedCell::no : { non_split_pixel(); break; } case SplittedCell::yes : { split_pixel(); break; } } }; auto constitutive_law_finite_strain = [this] (Strains_t Strains, Stresses_t Stresses, const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); // return value contains a tuple of rvalue_refs to both stress // and tangent moduli // for the case that cells do not contain splitt cells auto non_split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables, &grad](){ auto stress_tgt = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); Stresses = MatTB::PK1_stress (std::move(grad), std::move(std::get<0>(stress_tgt)), std::move( std::get<1>(stress_tgt))); }; // for the case that cells contain splitt cells auto split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables, &grad](){ auto stress_tgt = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); auto && stress_tgt_contributions = MatTB::PK1_stress (std::move(grad), std::move(std::get<0>(stress_tgt)), std::move(std::get<1>(stress_tgt))); auto && stress = std::get<0>(Stresses); auto && tangent = std::get<1>(Stresses); stress += ratio * std::get<0>(stress_tgt_contributions); tangent += ratio * std::get<1>(stress_tgt_contributions); }; switch (is_cell_splitted){ case SplittedCell::no : { non_split_pixel(); break; } case SplittedCell::yes : { split_pixel(); break; } } }; 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"); static_assert(std::is_same(arglist))>>::value, "Type mismatch for ratio reference, expected a real number"); 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,const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); // return value contains a tuple of rvalue_refs to both stress and tangent moduli auto & sigma = std::get<0>(Stresses); //for the case that cell does not contain any splitted pixel auto non_split_pixel = [&strain, &this_mat,ratio, &sigma, &internal_variables](){ sigma = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress(std::move(strain), internals...);}, internal_variables); }; // for the case that cells contain splitt cells auto split_pixel = [&strain, &this_mat, &ratio, &sigma, &internal_variables](){ sigma += ratio * apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress(std::move(strain), internals...);}, internal_variables); }; switch (is_cell_splitted){ case SplittedCell::no : { non_split_pixel(); break; } case SplittedCell::yes : { split_pixel(); break; } } }; auto constitutive_law_finite_strain = [this] (Strains_t Strains, Stresses_t && Stresses,const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); // 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. //for the case that cell does not contain any splitted pixel auto non_split_pixel = [ &strain, &this_mat, &ratio, &Stresses, &internal_variables, &grad] (){ 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 (grad, stress); }; // for the case that cells contain splitted cells auto split_pixel = [&strain, &this_mat, ratio, &Stresses, internal_variables, &grad](){ 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 += ratio * MatTB::PK1_stress (grad, stress); }; switch (is_cell_splitted){ case SplittedCell::no : { non_split_pixel(); break; } case SplittedCell::yes : { split_pixel(); break; } } }; iterable_proxy fields{*this, F, P}; for (auto && arglist: fields) { /** * arglist is a tuple of three tuples containing only Lvalue * references (see value_type 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(arglist))>>::value, "Type mismatch for ratio reference, expected a real number"); 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; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ 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(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.get_internals()}{}; /** 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(MaterialMuSpectre & mat, const StrainField_t & F, std::enable_if_t & P) :material{mat}, strain_field{F}, stress_tup{P}, internals{material.get_internals()}{}; //! Expected type for strain fields using StrainMap_t = typename traits::StrainMap_t; //! Expected type for stress fields using StressMap_t = typename traits::StressMap_t; //! Expected type for tangent stiffness fields using TangentMap_t = typename traits::TangentMap_t; //! expected type for strain values using Strain_t = typename traits::StrainMap_t::reference; //! expected type for stress values using Stress_t = typename traits::StressMap_t::reference; //! expected type for tangent stiffness values using Tangent_t = typename traits::TangentMap_t::reference; //! tuple of intervnal variables, depends on the material using InternalVariables = typename traits::InternalVariables; //! tuple containing a stress and possibly a tangent stiffness field using StressFieldTup = std::conditional_t <(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! tuple containing a stress and possibly a tangent stiffness field map using StressMapTup = std::conditional_t <(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! tuple containing a stress and possibly a tangent stiffness value ref 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) = default; //! Destructor virtual ~iterable_proxy() = default; //! Copy assignment operator iterable_proxy& operator=(const iterable_proxy &other) = default; //! Move assignment operator iterable_proxy& operator=(iterable_proxy &&other) = default; /** * dereferences into a tuple containing strains, and internal * variables, as well as maps to the stress and potentially * stiffness maps where to write the response of a pixel */ class iterator { public: //! type to refer to internal variables owned by a CRTP material using InternalReferences = MatTB::ReferenceTuple_t; //! return type to be unpacked per pixel my the constitutive law using value_type = std::tuple, Stress_tTup, Real, InternalReferences>; using iterator_category = std::forward_iterator_tag; //!< stl conformance //! 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) = default; //! Destructor virtual ~iterator() = 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; //!< ref to the proxy StrainMap_t strain_map; //!< map onto the global strain field //! map onto the global stress field and possibly tangent stiffness StressMapTup stress_map; size_t index; //!< index or pixel currently referred to private: }; //! returns iterator to first pixel if this material iterator begin() {return std::move(iterator(*this));} //! returns iterator past the last pixel in this material iterator end() {return std::move(iterator(*this, false));} protected: MaterialMuSpectre & material; //!< reference to the proxied material const StrainField_t & strain_field; //!< cell's global strain field //! references to the global stress field and perhaps tangent stiffness StressFieldTup stress_tup; //! references to the internal variables 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 && ratio = this->it.material.get_assigned_ratio(pixel); auto && stresses = apply([&pixel] (auto && ... stress_tgt) { return std::make_tuple(stress_tgt[pixel]...);}, this->stress_map); auto && internal = this->it.material.get_internals(); const auto id{this->index}; auto && internals = apply([id] (auto && ... internals_) { return InternalReferences{internals_[id]...};}, internal); return std::make_tuple(std::move(strain), std::move(stresses), std::move(ratio), std::move(internals)); } } // muSpectre #endif /* MATERIAL_MUSPECTRE_BASE_H */ diff --git a/src/materials/material_orthotropic.cc b/src/materials/material_orthotropic.cc new file mode 100644 index 0000000..46c1357 --- /dev/null +++ b/src/materials/material_orthotropic.cc @@ -0,0 +1,110 @@ +/** + * @file material_anisotropic.cc + * + * @author Ali Falsafi + * + * @date 11 Jul 2018 + * + * @brief Base class for materials (constitutive models) + * + * 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 "material_base.hh" +#include "common/common.hh" +#include "material_anisotropic.hh" +#include "material_orthotropic.hh" +#include "common/T4_map_proxy.hh" +//#include "common/tensor_algebra.hh" +//#include "common/eigen_tools.hh" + + +namespace muSpectre { + /* ---------------------------------------------------------------------- */ + template + MaterialOrthotropic::MaterialOrthotropic(std::string name, + std::vector input) + :Parent(name, input_c_maker(input)) + {}; + + template + std::vector + MaterialOrthotropic::input_c_maker(std::vector input){ + std::vector retval{}; + if (DimS==2){ + try{ + if (input.size() != 4){ + throw 1; + } + } + catch(int x) { + std::cout<<"Number of the entries should be 4"<; + //template class MaterialOrthotropic; + template class MaterialOrthotropic; +}//muSpectre diff --git a/src/materials/material_orthotropic.hh b/src/materials/material_orthotropic.hh new file mode 100644 index 0000000..96cbe8f --- /dev/null +++ b/src/materials/material_orthotropic.hh @@ -0,0 +1,135 @@ +/** + * @file material_anisotropic.hh + * + * @author Ali Falsafi + * + * @date 11 Jul 2018 + * + * @brief Base class for materials (constitutive models) + * + * 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_ORTHOTROPIC_H +#define MATERIAL_ORTHOTROPIC_H + +#include "materials/material_base.hh" +#include "materials/material_muSpectre_base.hh" +#include "materials/material_anisotropic.hh" +#include "common/common.hh" +#include "common/T4_map_proxy.hh" +#include "cell/cell_base.hh" + +namespace muSpectre { + + + // Forward declaration for factory function + //template + //class CellBase; + + template + class MaterialOrthotropic; + + //traits for orthotropic material + template + struct MaterialMuSpectre_traits> { + using Parent = MaterialMuSpectre_traits;//!< base for elasticity + + //! global field collection + using GFieldCollection_t = typename + MaterialBase::GFieldCollection_t; + + //! expected map type for strain fields + using StrainMap_t = MatrixFieldMap; + //! expected map type for stress fields + using StressMap_t = MatrixFieldMap; + //! expected map type for tangent stiffness fields + using TangentMap_t = T4MatrixFieldMap; + + //! declare what type of strain measure your law takes as input + constexpr static auto strain_measure{StrainMeasure::GreenLagrange}; + //! declare what type of stress measure your law yields as output + constexpr static auto stress_measure{StressMeasure::PK2}; + + //! anisotropicity without internal variables + using InternalVariables = std::tuple<>; + }; + /** + * Material implementation for orthotropic constitutive law + */ + template + class MaterialOrthotropic: + public MaterialAnisotropic + { + + //! base class + using Parent = MaterialAnisotropic; + + using Stiffness_t =T4Mat; + + //! traits of this material + using traits = MaterialMuSpectre_traits; + + //! this law does not have any internal variables + using InternalVariables = typename traits::InternalVariables; + + public: + //! global field collection + using GFieldCollection_t = typename + MaterialBase::GFieldCollection_t; + + //! expected map type for tangent stiffness fields + using Tangent_t = T4MatrixFieldMap; + + //! Default constructor + MaterialOrthotropic() = delete; + + MaterialOrthotropic(std::string name, std::vector input); + + //! Copy constructor + MaterialOrthotropic(const MaterialOrthotropic &other) = delete; + + //! Move constructor + MaterialOrthotropic(MaterialOrthotropic &&other) = delete; + + //! Destructor + virtual ~MaterialOrthotropic() = default; + + /* overloaded make function in order to make pyrhon binding + able to make an object of materila orthotropic*/ + static MaterialOrthotropic & make(CellBase & cell, + std::string name, std::vector input); + + std::vector input_c_maker(std::vector input); + + protected: + }; + + /* ---------------------------------------------------------------------- */ + template + MaterialOrthotropic & + MaterialOrthotropic::make(CellBase & cell, + std::string name, + std::vector input){ + auto mat = std::make_unique>(name, input); + auto & mat_ref = *mat; + cell.add_material(std::move(mat)); + return mat_ref; + } + +}//muSpectre +#endif diff --git a/src/materials/materials_toolbox.hh b/src/materials/materials_toolbox.hh index 3e2bd83..3243fd1 100644 --- a/src/materials/materials_toolbox.hh +++ b/src/materials/materials_toolbox.hh @@ -1,680 +1,693 @@ /** * @file materials_toolbox.hh * * @author Till Junge * * @date 02 Nov 2017 * * @brief collection of common continuum mechanics tools * * 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 MATERIALS_TOOLBOX_H #define MATERIALS_TOOLBOX_H #include "common/common.hh" #include "common/tensor_algebra.hh" #include "common/eigen_tools.hh" #include "common/T4_map_proxy.hh" #include #include #include #include #include #include #include namespace muSpectre { namespace MatTB { /** * thrown when generic materials-related runtime errors occur * (mostly continuum mechanics problems) */ class MaterialsToolboxError:public std::runtime_error{ public: //! constructor explicit MaterialsToolboxError(const std::string& what) :std::runtime_error(what){} //! constructor explicit MaterialsToolboxError(const char * what) :std::runtime_error(what){} }; /* ---------------------------------------------------------------------- */ /** * Flag used to designate whether the material should compute both stress * and tangent moduli or only stress */ enum class NeedTangent { yes, //!< compute both stress and tangent moduli no //!< compute only stress }; /** * struct used to determine the exact type of a tuple of references obtained * when a bunch of iterators over fiel_maps are dereferenced and their * results are concatenated into a tuple */ template struct ReferenceTuple { //! use this type using type = std::tuple; }; /** * specialisation for tuples */ //template <> template struct ReferenceTuple> { //! use this type using type = typename ReferenceTuple::type; }; /** * helper type for ReferenceTuple */ template using ReferenceTuple_t = typename ReferenceTuple::type; /* ---------------------------------------------------------------------- */ namespace internal { /** Structure for functions returning one strain measure as a * function of another **/ template struct ConvertStrain { static_assert((In == StrainMeasure::Gradient) || (In == StrainMeasure::Infinitesimal), "This situation makes me suspect that you are not using " "MatTb as intended. Disable this assert only if you are " "sure about what you are doing."); //! returns the converted strain template inline static decltype(auto) compute(Strain_t&& input) { // transparent case, in which no conversion is required: // just a perfect forwarding static_assert ((In == Out), "This particular strain conversion is not implemented"); return std::forward(input); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for getting Green-Lagrange strain from the transformation gradient E = ¹/₂ (C - I) = ¹/₂ (Fᵀ·F - I) **/ template <> struct ConvertStrain { //! returns the converted strain template inline static decltype(auto) compute(Strain_t&& F) { return .5*(F.transpose()*F - Strain_t::PlainObject::Identity()); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for getting Left Cauchy-Green strain from the transformation gradient B = F·Fᵀ = V² **/ template <> struct ConvertStrain { //! returns the converted strain template inline static decltype(auto) compute(Strain_t&& F) { return F*F.transpose(); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for getting Right Cauchy-Green strain from the transformation gradient C = Fᵀ·F = U² **/ template <> struct ConvertStrain { //! returns the converted strain template inline static decltype(auto) compute(Strain_t&& F) { return F.transpose()*F; } }; /* ---------------------------------------------------------------------- */ /** Specialisation for getting logarithmic (Hencky) strain from the transformation gradient E₀ = ¹/₂ ln C = ¹/₂ ln (Fᵀ·F) **/ template <> struct ConvertStrain { //! returns the converted strain template inline static decltype(auto) compute(Strain_t&& F) { constexpr Dim_t dim{EigenCheck::tensor_dim::value}; return (.5*logm(Eigen::Matrix{F.transpose()*F})).eval(); } }; } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning one strain measure as a function of //! another template decltype(auto) convert_strain(Strain_t && strain) { return internal::ConvertStrain::compute(std::move(strain)); }; /* ---------------------------------------------------------------------- */ namespace internal { /** Structure for functions returning PK1 stress from other stress measures **/ template struct PK1_stress { //! returns the converted stress template inline static decltype(auto) compute(Strain_t && /*strain*/, Stress_t && /*stress*/) { // the following test always fails to generate a compile-time error static_assert((StressM == StressMeasure::Cauchy) && (StressM == StressMeasure::PK1), "The requested Stress conversion is not implemented. " "You either made a programming mistake or need to " "implement it as a specialisation of this function. " "See PK2stress for an example."); } //! returns the converted stress and stiffness template inline static decltype(auto) compute(Strain_t && /*strain*/, Stress_t && /*stress*/, Tangent_t && /*stiffness*/) { // the following test always fails to generate a compile-time error static_assert((StressM == StressMeasure::Cauchy) && (StressM == StressMeasure::PK1), "The requested Stress conversion is not implemented. " "You either made a programming mistake or need to " "implement it as a specialisation of this function. " "See PK2stress for an example."); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for the transparent case, where we already have PK1 stress **/ template struct PK1_stress: public PK1_stress { //! returns the converted stress template inline static decltype(auto) compute(Strain_t && /*dummy*/, Stress_t && P) { return std::forward(P); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for the transparent case, where we already have PK1 stress *and* stiffness is given with respect to the transformation gradient **/ template struct PK1_stress: public PK1_stress { //! base class using Parent = PK1_stress; using Parent::compute; //! returns the converted stress and stiffness template inline static decltype(auto) compute(Strain_t && /*dummy*/, Stress_t && P, Tangent_t && K) { return std::make_tuple(std::forward(P), std::forward(K)); } }; /* ---------------------------------------------------------------------- */ /** * Specialisation for the case where we get material stress (PK2) */ template struct PK1_stress: public PK1_stress { //! returns the converted stress template inline static decltype(auto) compute(Strain_t && F, Stress_t && S) { return F*S; } }; /* ---------------------------------------------------------------------- */ /** * Specialisation for the case where we get material stress (PK2) derived * with respect to Green-Lagrange strain */ template struct PK1_stress: public PK1_stress { //! base class using Parent = PK1_stress; using Parent::compute; //! returns the converted stress and stiffness template inline static decltype(auto) compute(Strain_t && F, Stress_t && S, Tangent_t && C) { using T4 = typename std::remove_reference_t::PlainObject; using Tmap = T4MatMap; T4 K; Tmap Kmap{K.data()}; K.setZero(); for (int i = 0; i < Dim; ++i) { for (int m = 0; m < Dim; ++m) { for (int n = 0; n < Dim; ++n) { get(Kmap,i,m,i,n) += S(m,n); for (int j = 0; j < Dim; ++j) { for (int r = 0; r < Dim; ++r) { for (int s = 0; s < Dim; ++s) { get(Kmap,i,m,j,n) += F(i,r)*get(C,r,m,n,s)*(F(j,s)); } } } } } } auto && P = compute(std::forward(F), std::forward(S)); return std::make_tuple(std::move(P), std::move(K)); } }; /* ---------------------------------------------------------------------- */ /** * Specialisation for the case where we get Kirchhoff stress (τ) */ template struct PK1_stress: public PK1_stress { //! returns the converted stress template inline static decltype(auto) compute(Strain_t && F, Stress_t && tau) { return tau*F.inverse().transpose(); } }; /* ---------------------------------------------------------------------- */ /** * Specialisation for the case where we get Kirchhoff stress (τ) derived * with respect to Gradient */ template struct PK1_stress: public PK1_stress { //! short-hand using Parent = PK1_stress; using Parent::compute; //! returns the converted stress and stiffness template inline static decltype(auto) compute(Strain_t && F, Stress_t && tau, Tangent_t && C) { using T4 = typename std::remove_reference_t::PlainObject; using Tmap = T4MatMap; T4 K; Tmap Kmap{K.data()}; K.setZero(); auto && F_inv{F.inverse()}; for (int i = 0; i < Dim; ++i) { for (int m = 0; m < Dim; ++m) { for (int n = 0; n < Dim; ++n) { for (int j = 0; j < Dim; ++j) { for (int r = 0; r < Dim; ++r) { for (int s = 0; s < Dim; ++s) { get(Kmap,i,m,j,n) += F_inv(i,r)*get(C,r,m,n,s); } } } } } } auto && P = tau * F_inv.transpose(); return std::make_tuple(std::move(P), std::move(K)); } }; } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning an expression for PK2 stress based on template decltype(auto) PK1_stress(Strain_t && strain, Stress_t && stress) { constexpr Dim_t dim{EigenCheck::tensor_dim::value}; static_assert((dim == EigenCheck::tensor_dim::value), "Stress and strain tensors have differing dimensions"); return internal::PK1_stress::compute (std::forward(strain), std::forward(stress)); }; /* ---------------------------------------------------------------------- */ //! set of functions returning an expression for PK2 stress based on template decltype(auto) PK1_stress(Strain_t && strain, Stress_t && stress, Tangent_t && tangent) { constexpr Dim_t dim{EigenCheck::tensor_dim::value}; static_assert((dim == EigenCheck::tensor_dim::value), "Stress and strain tensors have differing dimensions"); static_assert((dim == EigenCheck::tensor_4_dim::value), "Stress and tangent tensors have differing dimensions"); return internal::PK1_stress::compute (std::forward(strain), std::forward(stress), std::forward(tangent)); }; namespace internal { //! Base template for elastic modulus conversion template struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real& /*in1*/, const Real& /*in2*/) { // if no return has happened until now, the conversion is not // implemented yet static_assert((In1 == In2), "This conversion has not been implemented yet, please add " "it here below as a specialisation of this function " "template. Check " "https://en.wikipedia.org/wiki/Lam%C3%A9_parameters for " "// TODO: he formula."); return 0; } }; /** * Spectialisation for when the output is the first input */ template struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & A, const Real & /*B*/) { return A; } }; /** * Spectialisation for when the output is the second input */ template struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & /*A*/, const Real & B) { return B; } }; /** * Specialisation μ(E, ν) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & E, const Real & nu) { return E/(2*(1+nu)); } }; /** * Specialisation λ(E, ν) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & E, const Real & nu) { return E*nu/((1+nu)*(1-2*nu)); } }; /** * Specialisation K(E, ν) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & E, const Real & nu) { return E/(3*(1-2*nu)); } }; /** * Specialisation E(K, µ) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & K, const Real & G) { return 9*K*G/(3*K+G); } }; /** * Specialisation ν(K, µ) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & K, const Real & G) { return (3*K - 2*G) /(2*(3*K + G)); } }; /** * Specialisation E(λ, µ) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & lambda, const Real & G) { return G * (3*lambda + 2*G)/(lambda + G); } }; } // internal /** * allows the conversion from any two distinct input moduli to a * chosen output modulus */ template inline constexpr Real convert_elastic_modulus(const Real& in1, const Real& in2) { // enforcing sanity static_assert((In1 != In2), "The input modulus types cannot be identical"); // enforcing independence from order in which moduli are supplied constexpr bool inverted{In1 > In2}; using Converter = std::conditional_t, internal::Converter>; if (inverted) { return Converter::compute(std::move(in2), std::move(in1)); } else { return Converter::compute(std::move(in1), std::move(in2)); } } //! static inline implementation of Hooke's law template struct Hooke { /** * compute Lamé's first constant * @param young: Young's modulus * @param poisson: Poisson's ratio */ inline static constexpr Real compute_lambda(const Real & young, const Real & poisson) { return convert_elastic_modulus(young, poisson); } /** * compute Lamé's second constant (i.e., shear modulus) * @param young: Young's modulus * @param poisson: Poisson's ratio */ inline static constexpr Real compute_mu(const Real & young, const Real & poisson) { return convert_elastic_modulus(young, poisson); } /** * compute the bulk modulus * @param young: Young's modulus * @param poisson: Poisson's ratio */ inline static constexpr Real compute_K(const Real & young, const Real & poisson) { return convert_elastic_modulus(young, poisson); } /** * compute the stiffness tensor * @param lambda: Lamé's first constant * @param mu: Lamé's second constant (i.e., shear modulus) */ inline static Eigen::TensorFixedSize> compute_C(const Real & lambda, const Real & mu) { return lambda*Tensors::outer(Tensors::I2(),Tensors::I2()) + 2*mu*Tensors::I4S(); } /** * compute the stiffness tensor * @param lambda: Lamé's first constant * @param mu: Lamé's second constant (i.e., shear modulus) */ inline static T4Mat compute_C_T4(const Real & lambda, const Real & mu) { return lambda*Matrices::Itrac() + 2*mu*Matrices::Isymm(); } /** * return stress * @param lambda: First Lamé's constant * @param mu: Second Lamé's constant (i.e. shear modulus) * @param E: Green-Lagrange or small strain tensor */ template inline static decltype(auto) evaluate_stress(const Real & lambda, const Real & mu, s_t && E) { return E.trace()*lambda * Strain_t::Identity() + 2*mu*E; } + + /** + * return stress + * @param C: stiffness tensor (Piola-Kirchhoff 2 (or σ) w.r.t to `E`) + * @param E: Green-Lagrange or small strain tensor + */ + /*template + inline static decltype(auto) + evaluate_stress(const Tangent_t C, s_t && E) { + return tensmult(C, E); + }*/ + + /** * return stress and tangent stiffness * @param lambda: First Lamé's constant * @param mu: Second Lamé's constant (i.e. shear modulus) * @param E: Green-Lagrange or small strain tensor * @param C: stiffness tensor (Piola-Kirchhoff 2 (or σ) w.r.t to `E`) */ template inline static decltype(auto) evaluate_stress(const Real & lambda, const Real & mu, Tangent_t && C, s_t && E) { return std::make_tuple (std::move(evaluate_stress(lambda, mu, std::move(E))), std::move(C)); } }; } // MatTB } // muSpectre #endif /* MATERIALS_TOOLBOX_H */ diff --git a/tests/test_solver_newton_cg.cc b/tests/test_solver_newton_cg.cc index 88e31ed..62bcb60 100644 --- a/tests/test_solver_newton_cg.cc +++ b/tests/test_solver_newton_cg.cc @@ -1,447 +1,456 @@ /** * @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 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 "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 "fft/fftw_engine.hh" #include "fft/projection_finite_strain_fast.hh" #include "materials/material_linear_elastic1.hh" +#include "materials/material_orthotropic.hh" #include "common/iterators.hh" #include "common/ccoord_operations.hh" #include "common/common.hh" #include "cell/cell_factory.hh" #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 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}; + constexpr Ccoord_t resolutions{5, 5}; + constexpr Rcoord_t lengths{5, 5}; auto fft_ptr{std::make_unique>(resolutions, ipow(dim, 2))}; - auto proj_ptr{std::make_unique>(std::move(fft_ptr), lengths)}; + auto proj_ptr{std::make_unique>(std::move(fft_ptr), lengths)}; CellBase sys(std::move(proj_ptr)); - using Mat_t = MaterialLinearElastic1; + using Mat_t = MaterialOrthotropic; + //using Mat_t = MaterialLinearElastic1; //const Real Young{210e9}, Poisson{.33}; - const Real Young{1.0030648180242636}, Poisson{0.29930675909878679}; + const Real E{210.0e9}, n{0.3}; + const int con {10} ; // 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); + //auto& Material_hard = Mat_t::make(sys, "hard", 10*Young, Poisson); + //auto& Material_soft = Mat_t::make(sys, "soft", Young, Poisson); + + std::vector input_soft;input_soft.push_back((E*(1-n))/((1+n)*(1-2*n))); input_soft.push_back((E*(n))/((1+n)*(1-2*n))); input_soft.push_back((E*(1-n))/((1+n)*(1-2*n))); input_soft.push_back((E*(1-n))/((1+n)*(2)));//input_soft.push_back(0.0); input_soft.push_back(0.8); + std::vector input_hard;input_hard.push_back(con*(E*(1-n))/((1+n)*(1-2*n))); input_hard.push_back(con*(E*(n))/((1+n)*(1-2*n))); input_hard.push_back(con*(E*(1-n))/((1+n)*(1-2*n))); input_hard.push_back(con*(E)/((1+n)*(2)));//input_hard.push_back(0.0); input_hard.push_back(0.8); + + auto& Material_hard = Mat_t::make(sys, "hard", input_hard); + auto& Material_soft = Mat_t::make(sys, "soft", input_soft); for (auto && tup: akantu::enumerate(sys)) { auto && pixel = std::get<1>(tup); if (std::get<0>(tup) == 0) { Material_hard.add_pixel(pixel); - } else { + } else { Material_soft.add_pixel(pixel); - } + } } sys.initialise(); Grad_t delF0; - delF0 << 0, 1., 0, 0, 0, 0, 0, 0, 0; + delF0 << 1e-4, 0, 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}; GradIncrements grads; grads.push_back(delF0); DeprecatedSolverCG cg{sys, cg_tol, maxiter, bool(verbose)}; Eigen::ArrayXXd res1{deprecated_de_geus(sys, grads, cg, newton_tol, verbose)[0].grad}; DeprecatedSolverCG cg2{sys, cg_tol, maxiter, bool(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 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 Uint maxiter{dim*1000}; constexpr Dim_t verbose{0}; DeprecatedSolverCGEigen cg{sys, cg_tol, maxiter, bool(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, bool(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>; 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 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, bool(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, bool(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 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, bool(verbose)}; auto F = sys.get_strain_vector(); F.setZero(); sys.evaluate_stress_tangent(); Eigen::VectorXd DelF(sys.get_nb_dof()); using RMap_t = 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, bool(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(); } // muSpectre