diff --git a/language_bindings/python/bind_py_cell.cc b/language_bindings/python/bind_py_cell.cc index f40b79c..d428675 100644 --- a/language_bindings/python/bind_py_cell.cc +++ b/language_bindings/python/bind_py_cell.cc @@ -1,309 +1,307 @@ /** * @file bind_py_cell.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief Python bindings for the cell factory function * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/common.hh" #include "common/ccoord_operations.hh" #include "cell/cell_factory.hh" #include "cell/cell_split_factory.hh" #include "cell/cell_base.hh" #include "cell/cell_split.hh" #ifdef WITH_FFTWMPI #include "fft/fftwmpi_engine.hh" #endif #ifdef WITH_PFFT #include "fft/pfft_engine.hh" #endif #include #include #include #include "pybind11/eigen.h" #include #include using muSpectre::Ccoord_t; using muSpectre::Dim_t; using muSpectre::Formulation; using muSpectre::Rcoord_t; using pybind11::literals::operator""_a; namespace py = pybind11; /** * cell factory for specific FFT engine */ #ifdef WITH_MPI template void add_parallel_cell_factory_helper(py::module & mod, const char * name) { using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; mod.def(name, [](Ccoord res, Rcoord lens, Formulation form, size_t comm) { return make_parallel_cell, FFTEngine>( std::move(res), std::move(lens), std::move(form), std::move(Communicator(MPI_Comm(comm)))); }, "resolutions"_a, "lengths"_a = CcoordOps::get_cube(1.), "formulation"_a = Formulation::finite_strain, "communicator"_a = size_t(MPI_COMM_SELF)); } #endif /** * the cell factory is only bound for default template parameters */ template void add_cell_factory_helper(py::module & mod) { using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; mod.def("CellFactory", [](Ccoord res, Rcoord lens, Formulation form) { return muSpectre::make_cell_base_ptr_cell_base< dim, dim, muSpectre::CellBase>( std::move(res), std::move(lens), std::move(form)); }, "resolutions"_a, "lengths"_a = muSpectre::CcoordOps::get_cube(1.), "formulation"_a = Formulation::finite_strain); mod.def("CellFactorySplit", [](Ccoord res, Rcoord lens, Formulation form) { - // return make_cell_ptr>(std::move(res), std::move(lens), std::move(form)); return muSpectre::make_cell_base_ptr_cell_split< dim, dim, muSpectre::CellBase>( std::move(res), std::move(lens), std::move(form)); }, "resolutions"_a, "lengths"_a = muSpectre::CcoordOps::get_cube(1.), "formulation"_a = Formulation::finite_strain); #ifdef WITH_FFTWMPI add_parallel_cell_factory_helper>( mod, "FFTWMPICellFactory"); #endif #ifdef WITH_PFFT add_parallel_cell_factory_helper>(mod, "PFFTCellFactory"); #endif } void add_cell_factory(py::module & mod) { add_cell_factory_helper(mod); add_cell_factory_helper(mod); } /** * CellBase for which the material and spatial dimension are identical */ template void add_cell_base_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "CellBase" << dim << 'd'; const std::string name = name_stream.str(); using sys_t = muSpectre::CellBase; py::class_>(mod, name.c_str()) .def("__len__", &sys_t::size) .def("__iter__", [](sys_t & s) { return py::make_iterator(s.begin(), s.end()); }) .def("initialise", &sys_t::initialise, "flags"_a = muSpectre::FFT_PlanFlags::estimate) .def( "directional_stiffness", [](sys_t & cell, py::EigenDRef & v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim * dim)) { std::stringstream err{}; err << "need array of shape (" << dim * dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } if (!cell.is_initialised()) { cell.initialise(); } const std::string out_name{"temp output for directional stiffness"}; const std::string in_name{"temp input for directional stiffness"}; constexpr bool create_tangent{true}; auto & K = cell.get_tangent(create_tangent); auto & input = cell.get_managed_T2_field(in_name); auto & output = cell.get_managed_T2_field(out_name); input.eigen() = v; cell.directional_stiffness(K, input, output); return output.eigen(); }, "δF"_a) .def("project", [](sys_t & cell, py::EigenDRef & v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim * dim)) { std::stringstream err{}; err << "need array of shape (" << dim * dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } if (!cell.is_initialised()) { cell.initialise(); } const std::string in_name{"temp input for projection"}; auto & input = cell.get_managed_T2_field(in_name); input.eigen() = v; cell.project(input); return input.eigen(); }, "field"_a) .def("get_strain", [](sys_t & s) { return s.get_strain().eigen(); }, py::return_value_policy::reference_internal) .def("get_stress", [](sys_t & s) { return Eigen::ArrayXXd(s.get_stress().eigen()); }) .def_property_readonly("size", &sys_t::size) .def("evaluate_stress_tangent", [](sys_t & cell, py::EigenDRef & v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim * dim)) { std::stringstream err{}; err << "need array of shape (" << dim * dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } auto & strain{cell.get_strain()}; strain.eigen() = v; auto stress_tgt{cell.evaluate_stress_tangent(strain)}; return std::tuple( std::get<0>(stress_tgt).eigen(), std::get<1>(stress_tgt).eigen()); }, "strain"_a) .def("evaluate_stress", [](sys_t & cell, py::EigenDRef & v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim * dim)) { std::stringstream err{}; err << "need array of shape (" << dim * dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } auto & strain{cell.get_strain()}; strain.eigen() = v; return cell.evaluate_stress(); }, "strain"_a, py::return_value_policy::reference_internal) .def("get_projection", &sys_t::get_projection) .def("get_subdomain_resolutions", &sys_t::get_subdomain_resolutions) .def("get_subdomain_locations", &sys_t::get_subdomain_locations) .def("get_domain_resolutions", &sys_t::get_domain_resolutions) .def("get_domain_lengths", &sys_t::get_domain_resolutions) .def("set_uniform_strain", [](sys_t & cell, py::EigenDRef & v) -> void { cell.set_uniform_strain(v); }, "strain"_a) .def("save_history_variables", &sys_t::save_history_variables); } template void add_cell_split_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "CellSplit" << dim << 'd'; const auto name{name_stream.str()}; using sys_split_t = muSpectre::CellSplit; using sys_t = muSpectre::CellBase; using Mat_t = muSpectre::MaterialBase; py::class_>(mod, name.c_str()) .def("make_precipitate", [](sys_split_t & cell, Mat_t & mat, std::vector> precipitate_vertices) { cell.make_automatic_precipitate(precipitate_vertices, mat); }, "material"_a, "vertices"_a) .def("complete_material_assignment", [](sys_split_t & cell, Mat_t & mat) { cell.complete_material_assignment(mat); }, "material"_a); } void add_cell_base(py::module & mod) { py::class_>(mod, "Cell") .def("get_globalised_internal_real_array", &muSpectre::Cell::get_globalised_internal_real_array, "unique_name"_a, "Convenience function to copy local (internal) fields of " "materials into a global field. At least one of the materials in " "the cell needs to contain an internal field named " "`unique_name`. If multiple materials contain such a field, they " "all need to be of same scalar type and same number of " "components. This does not work for split pixel cells or " "laminate pixel cells, as they can have multiple entries for the " "same pixel. Pixels for which no field named `unique_name` " "exists get an array of zeros." "\n" "Parameters:\n" "unique_name: fieldname to fill the global field with. At " "least one material must have such a field, or an " "Exception is raised.") .def("get_globalised_current_real_array", &muSpectre::Cell::get_globalised_current_real_array, "unique_name"_a) .def("get_globalised_old_real_array", &muSpectre::Cell::get_globalised_old_real_array, "unique_name"_a, "nb_steps_ago"_a = 1) .def("get_managed_real_array", &muSpectre::Cell::get_managed_real_array, "unique_name"_a, "nb_components"_a, "returns a field or nb_components real numbers per pixel"); add_cell_base_helper(mod); add_cell_base_helper(mod); } void add_cell_split(py::module & mod) { add_cell_split_helper(mod); add_cell_split_helper(mod); } void add_cell(py::module & mod) { add_cell_factory(mod); auto cell{mod.def_submodule("cell")}; cell.doc() = "bindings for cells and cell factories"; cell.def("scale_by_2", [](py::EigenDRef & v) { v *= 2; }); add_cell_base(cell); add_cell_split(cell); } diff --git a/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc index 5f64d0b..9ab1b64 100644 --- a/language_bindings/python/bind_py_material.cc +++ b/language_bindings/python/bind_py_material.cc @@ -1,261 +1,269 @@ /** * @file bind_py_material.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief python bindings for µSpectre's materials * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/common.hh" #include "materials/material_anisotropic.hh" #include "materials/material_orthotropic.hh" #include "cell/cell_base.hh" #include "cell/cell_split.hh" #include "materials/material_base.hh" #include "materials/material_evaluator.hh" #include #include #include #include #include using muSpectre::Dim_t; using muSpectre::Real; using pybind11::literals::operator""_a; namespace py = pybind11; /* ---------------------------------------------------------------------- */ template void add_material_linear_elastic_generic1_helper(py::module & mod); /* ---------------------------------------------------------------------- */ template void add_material_linear_elastic_generic2_helper(py::module & mod); 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 = muSpectre::MaterialBase; py::class_(mod, name.c_str()); } /** * python binding for the optionally objective form of Hooke's law */ template void add_material_linear_elastic1_helper(py::module & mod); template void add_material_linear_elastic2_helper(py::module & mod); template void add_material_linear_elastic3_helper(py::module & mod); template void add_material_linear_elastic4_helper(py::module & mod); template void add_material_hyper_elasto_plastic1_helper(py::module & mod); template void add_material_anisotropic_helper(py::module & mod); template void add_material_orthotropic_helper(py::module & mod); /* ---------------------------------------------------------------------- */ template class PyMaterialBase : public muSpectre::MaterialBase { public: /* Inherit the constructors */ using Parent = muSpectre::MaterialBase; using Parent::Parent; /* Trampoline (need one for each virtual function) */ void save_history_variables() override { PYBIND11_OVERLOAD_PURE(void, // Return type Parent, // Parent class save_history_variables); // Name of function in C++ // (must match Python name) } /* Trampoline (need one for each virtual function) */ void initialise() override { PYBIND11_OVERLOAD_PURE( void, // Return type Parent, // Parent class initialise); // Name of function in C++ (must match Python name) } void compute_stresses(const typename Parent::StrainField_t & F, typename Parent::StressField_t & P, muSpectre::Formulation form, muSpectre::SplitCell is_cell_split) override { PYBIND11_OVERLOAD_PURE( void, // Return type Parent, // Parent class compute_stresses, // Name of function in C++ (must match Python name) F, P, form, is_cell_split); } void compute_stresses_tangent(const typename Parent::StrainField_t & F, typename Parent::StressField_t & P, typename Parent::TangentField_t & K, muSpectre::Formulation form, muSpectre::SplitCell is_cell_split) override { PYBIND11_OVERLOAD_PURE( void, /* Return type */ Parent, /* Parent class */ compute_stresses, /* Name of function in C++ (must match Python name) */ F, P, K, form, is_cell_split); } }; template void add_material_evaluator(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialEvaluator_" << Dim << "d"; std::string name{name_stream.str()}; using MatEval_t = muSpectre::MaterialEvaluator; py::class_(mod, name.c_str()) .def(py::init>>()) .def("save_history_variables", &MatEval_t::save_history_variables, "for materials with state variables") .def("evaluate_stress", [](MatEval_t & mateval, py::EigenDRef & grad, muSpectre::Formulation form) { if ((grad.cols() != Dim) or (grad.rows() != Dim)) { std::stringstream err{}; err << "need matrix of shape (" << Dim << "×" << Dim << ") but got (" << grad.rows() << "×" << grad.cols() << ")."; throw std::runtime_error(err.str()); } return mateval.evaluate_stress(grad, form); }, "strain"_a, "formulation"_a, "Evaluates stress for a given strain and formulation " "(Piola-Kirchhoff 1 stress as a function of the placement gradient " "P = P(F) for formulation=Formulation.finite_strain and Cauchy " "stress as a function of the infinitesimal strain tensor σ = σ(ε) " "for formulation=Formulation.small_strain)") .def("evaluate_stress_tangent", [](MatEval_t & mateval, py::EigenDRef & grad, muSpectre::Formulation form) { if ((grad.cols() != Dim) or (grad.rows() != Dim)) { std::stringstream err{}; err << "need matrix of shape (" << Dim << "×" << Dim << ") but got (" << grad.rows() << "×" << grad.cols() << ")."; throw std::runtime_error(err.str()); } return mateval.evaluate_stress_tangent(grad, form); }, "strain"_a, "formulation"_a, "Evaluates stress and tangent moduli for a given strain and " "formulation (Piola-Kirchhoff 1 stress as a function of the " "placement gradient P = P(F) for " "formulation=Formulation.finite_strain and Cauchy stress as a " "function of the infinitesimal strain tensor σ = σ(ε) for " "formulation=Formulation.small_strain). The tangent moduli are K = " "∂P/∂F for formulation=Formulation.finite_strain and C = ∂σ/∂ε for " "formulation=Formulation.small_strain.") .def("estimate_tangent", [](MatEval_t & evaluator, py::EigenDRef & grad, muSpectre::Formulation form, const Real step, const muSpectre::FiniteDiff diff_type) { if ((grad.cols() != Dim) or (grad.rows() != Dim)) { std::stringstream err{}; err << "need matrix of shape (" << Dim << "×" << Dim << ") but got (" << grad.rows() << "×" << grad.cols() << ")."; throw std::runtime_error(err.str()); } return evaluator.estimate_tangent(grad, form, step, diff_type); }, "strain"_a, "formulation"_a, "Delta_x"_a, "difference_type"_a = muSpectre::FiniteDiff::centred, "Numerical estimate of the tangent modulus using finite " "differences. The finite difference scheme as well as the finite " "step size can be chosen. If there are no special circumstances, " "the default scheme of centred finite differences yields the most " "accurate results at an increased computational cost."); } template void add_material_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialBase_" << Dim << "d"; std::string name{name_stream.str()}; using Material = muSpectre::MaterialBase; using MaterialTrampoline = PyMaterialBase; using FC_t = muSpectre::LocalFieldCollection; using FCBase_t = muSpectre::FieldCollectionBase; py::class_>(mod, name.c_str()) .def(py::init()) .def("save_history_variables", &Material::save_history_variables) .def("list_fields", &Material::list_fields) .def("get_real_field", &Material::get_real_field, "field_name"_a, py::return_value_policy::reference_internal) .def("size", &Material::size) .def("add_pixel", [](Material & mat, muSpectre::Ccoord_t pix) { mat.add_pixel(pix); }, "pixel"_a) .def("add_pixel_split", [](Material & mat, muSpectre::Ccoord_t pix, muSpectre::Real ratio) { mat.add_pixel_split(pix, ratio); }, "pixel"_a, "ratio"_a) + + .def_static( + "make_C_Voigt_from_col", + [](Eigen::Matrix & C_col) { + return Material::make_C_voigt_from_col(C_col); + }, + "C_col"_a) + .def_property_readonly("collection", [](Material & material) -> FCBase_t & { return material.get_collection(); }, "returns the field collection containing internal " "fields of this material"); add_material_linear_elastic1_helper(mod); add_material_linear_elastic2_helper(mod); add_material_linear_elastic3_helper(mod); add_material_linear_elastic4_helper(mod); add_material_hyper_elasto_plastic1_helper(mod); add_material_linear_elastic_generic1_helper(mod); add_material_linear_elastic_generic2_helper(mod); add_material_anisotropic_helper(mod); add_material_orthotropic_helper(mod); add_material_evaluator(mod); } void add_material(py::module & mod) { auto material{mod.def_submodule("material")}; material.doc() = "bindings for constitutive laws"; add_material_helper(material); add_material_helper(material); } diff --git a/language_bindings/python/bind_py_material_linear_elastic1.cc b/language_bindings/python/bind_py_material_linear_elastic1.cc index cfe31d0..383fa09 100644 --- a/language_bindings/python/bind_py_material_linear_elastic1.cc +++ b/language_bindings/python/bind_py_material_linear_elastic1.cc @@ -1,72 +1,73 @@ /** * @file bind_py_material_linear_elastic1.cc * * @author Till Junge * * @date 29 Oct 2018 * * @brief python bindings for MaterialLinearElastic1 * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "materials/material_linear_elastic1.hh" #include "cell/cell_base.hh" #include #include #include #include #include using muSpectre::Dim_t; using muSpectre::Real; using pybind11::literals::operator""_a; namespace py = pybind11; /** * 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 = muSpectre::MaterialLinearElastic1; using Sys_t = muSpectre::CellBase; py::class_, std::shared_ptr>( 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_static("make_evaluator", [](Real e, Real p) { return Mat_t::make_evaluator(e, p); }, - "Young"_a, "Poisson"_a); + "Young"_a, "Poisson"_a) + .def("get_C", [](Mat_t & mat) { return mat.get_C(); }); } template void add_material_linear_elastic1_helper(py::module &); template void add_material_linear_elastic1_helper(py::module &); diff --git a/language_bindings/python/bind_py_material_linear_elastic_generic.cc b/language_bindings/python/bind_py_material_linear_elastic_generic.cc index 9ca6771..8a5ebf0 100644 --- a/language_bindings/python/bind_py_material_linear_elastic_generic.cc +++ b/language_bindings/python/bind_py_material_linear_elastic_generic.cc @@ -1,150 +1,151 @@ /** * @file bind_py_material_linear_elastic_generic.cc * * @author Till Junge * * @date 20 Dec 2018 * * @brief bindings for the generic linear elastic law defined by its stiffness * tensor * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "materials/material_linear_elastic_generic1.hh" #include "materials/material_linear_elastic_generic2.hh" #include "cell/cell_base.hh" #include #include #include #include #include using namespace muSpectre; // NOLINT // TODO(junge): figure this out namespace py = pybind11; using namespace pybind11::literals; // NOLINT: recommended usage /** * python binding for the generic linear elastic material */ template void add_material_linear_elastic_generic1_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElasticGeneric1_" << Dim << "d"; const auto name{name_stream.str()}; using Mat_t = MaterialLinearElasticGeneric1; using Cell_t = CellBase; py::class_, std::shared_ptr>( mod, name.c_str()) .def_static( "make", [](Cell_t & cell, std::string name, const py::EigenDRef< Eigen::Matrix> & elastic_tensor) -> Mat_t & { return Mat_t::make(cell, name, elastic_tensor); }, "cell"_a, "name"_a, "elastic_tensor"_a, py::return_value_policy::reference, py::keep_alive<1, 0>(), "Factory function returning a MaterialLinearElastic instance. " "The elastic tensor has to be specified in Voigt notation.") .def("add_pixel", [](Mat_t & mat, Ccoord_t pix) { mat.add_pixel(pix); }, "pixel"_a, "Register a new pixel to this material. subsequent evaluations of " "the " "stress and tangent in the cell will use this constitutive law for " "this " "particular pixel") .def("size", &Mat_t::size) .def_static("make_evaluator", [](const py::EigenDRef< Eigen::Matrix> & elastic_tensor) { return Mat_t::make_evaluator(elastic_tensor); }, - "elastic_tensor"_a); + "elastic_tensor"_a) + .def("get_C", [](Mat_t & mat) { return mat.get_C(); }); } /** * python binding for the generic linear elastic material with eigenstcain */ template void add_material_linear_elastic_generic2_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElasticGeneric2_" << Dim << "d"; const auto name{name_stream.str()}; using Mat_t = MaterialLinearElasticGeneric2; using Cell_t = CellBase; py::class_, std::shared_ptr>( mod, name.c_str()) .def_static( "make", [](Cell_t & cell, std::string name, const py::EigenDRef< Eigen::Matrix> & elastic_tensor) -> Mat_t & { return Mat_t::make(cell, name, elastic_tensor); }, "cell"_a, "name"_a, "elastic_tensor"_a, py::return_value_policy::reference, py::keep_alive<1, 0>(), "Factory function returning a MaterialLinearElastic instance. " "The elastic tensor has to be specified in Voigt notation.") .def("add_pixel", [](Mat_t & mat, Ccoord_t pix) { mat.add_pixel(pix); }, "pixel"_a, "Register a new pixel to this material. Subsequent evaluations of " "the stress and tangent in the cell will use this constitutive law " "for this particular pixel") .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, "Register a new pixel to this material and assign the eigenstrain. " "Subsequent Evaluations of the stress and tangent in the cell will " "use this constitutive law for this particular pixel") .def("size", &Mat_t::size) .def_static("make_evaluator", [](const py::EigenDRef< Eigen::Matrix> & elastic_tensor) { return Mat_t::make_evaluator(elastic_tensor); }, "elastic_tensor"_a); } template void add_material_linear_elastic_generic1_helper(py::module &); template void add_material_linear_elastic_generic1_helper(py::module &); template void add_material_linear_elastic_generic2_helper(py::module &); template void add_material_linear_elastic_generic2_helper(py::module &); diff --git a/src/materials/material_anisotropic.hh b/src/materials/material_anisotropic.hh index 8f9e05a..5beb366 100644 --- a/src/materials/material_anisotropic.hh +++ b/src/materials/material_anisotropic.hh @@ -1,156 +1,156 @@ /** * @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 SRC_MATERIALS_MATERIAL_ANISOTROPIC_HH_ #define SRC_MATERIALS_MATERIAL_ANISOTROPIC_HH_ #include "materials/material_base.hh" #include "materials/material_muSpectre_base.hh" #include "materials/stress_transformations_PK2.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, 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: EIGEN_MAKE_ALIGNED_OPERATOR_NEW 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; // constructor // a std::vector is utilized as the input of the constructor to // enable us to check its length so to prevent user mistake 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: static auto c_maker(std::vector input) -> Stiffness_t; - auto get_stiffness() { return this->C; } + auto get_C() { return this->C; } 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); } } // namespace muSpectre #endif // SRC_MATERIALS_MATERIAL_ANISOTROPIC_HH_ diff --git a/src/materials/material_base.cc b/src/materials/material_base.cc index 4c54711..069ed78 100644 --- a/src/materials/material_base.cc +++ b/src/materials/material_base.cc @@ -1,172 +1,191 @@ /** * @file material_base.cc * * @author Till Junge * * @date 01 Nov 2017 * * @brief implementation of material * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "materials/material_base.hh" namespace muSpectre { //----------------------------------------------------------------------------// template MaterialBase::MaterialBase(std::string name) : name(name) { static_assert((DimM == oneD) || (DimM == twoD) || (DimM == threeD), "only 1, 2, or threeD supported"); } /* ---------------------------------------------------------------------- */ template const std::string & MaterialBase::get_name() const { return this->name; } /* ---------------------------------------------------------------------- */ template void MaterialBase::add_pixel(const Ccoord & ccoord) { this->internal_fields.add_pixel(ccoord); } /* ---------------------------------------------------------------------- */ template void MaterialBase::add_pixel_split(const Ccoord & local_ccoord, Real ratio) { auto & this_mat = static_cast(*this); this_mat.add_pixel(local_ccoord); this->assigned_ratio.value().get().push_back(ratio); } /* ---------------------------------------------------------------------- */ template void MaterialBase::add_pixel_laminate(const Ccoord & local_ccoord, Real ratio, Rcoord normal_vec) { auto & this_mat = static_cast(*this); this_mat.add_pixel(local_ccoord); this->assigned_ratio.value().get().push_back(ratio); Eigen::Map> normal_vec_map(normal_vec.data()); this->assigned_normal_vec.value().get().push_back(normal_vec_map); } /* ---------------------------------------------------------------------- */ template void MaterialBase::allocate_optional_fields(SplitCell is_cell_split) { if (is_cell_split == SplitCell::simple || is_cell_split == SplitCell::laminate) { this->assigned_ratio = make_field("ratio", this->internal_fields); if (is_cell_split == SplitCell::laminate) { this->assigned_normal_vec = make_field("normal vector", this->internal_fields); // this->local_F = // make_field("local Strain", // this->internal_fields); // this->local_P = // make_field("local Stress", // this->internal_fields); // this->local_K = // make_field("local Tangent", // this->internal_fields); } } } /* ---------------------------------------------------------------------- */ template void MaterialBase::compute_stresses(const Field_t & F, Field_t & P, Formulation form, SplitCell is_cell_split) { this->compute_stresses(StrainField_t::check_ref(F), StressField_t::check_ref(P), form, is_cell_split); } /* ---------------------------------------------------------------------- */ template void MaterialBase::get_assigned_ratios( std::vector & pixel_assigned_ratios, Ccoord subdomain_resolutions, Ccoord subdomain_locations) { auto assigned_ratio_mapped = this->assigned_ratio.value().get().get_map(); for (auto && key_val : this->assigned_ratio.value().get().get_map().enumerate()) { const auto & ccoord = std::get<0>(key_val); const auto & val = std::get<1>(key_val); auto index = CcoordOps::get_index(subdomain_resolutions, subdomain_locations, ccoord); pixel_assigned_ratios[index] += val; } } /* ---------------------------------------------------------------------- */ template Real MaterialBase::get_assigned_ratio(Ccoord pixel) { return this->assigned_ratio.value().get().get_map()[pixel]; } //----------------------------------------------------------------------------// template void MaterialBase::compute_stresses_tangent( const Field_t & F, Field_t & P, Field_t & K, Formulation form, SplitCell is_cell_split) { this->compute_stresses_tangent( StrainField_t::check_ref(F), StressField_t::check_ref(P), TangentField_t::check_ref(K), form, is_cell_split); } template auto MaterialBase::get_real_field(std::string field_name) -> EigenMap { if (not this->internal_fields.check_field_exists(field_name)) { std::stringstream err{}; err << "Field '" << field_name << "' does not exist in material '" << this->name << "'."; throw FieldCollectionError(err.str()); } auto & field{this->internal_fields[field_name]}; if (field.get_stored_typeid().hash_code() != typeid(Real).hash_code()) { std::stringstream err{}; err << "Field '" << field_name << "' is not real-valued"; throw FieldCollectionError(err.str()); } return static_cast &>(field).eigen(); } + /* ---------------------------------------------------------------------- */ + template + auto MaterialBase::make_C_voigt_from_col( + T4Mat C_col) -> CVoigt_t { + using VC_t = VoigtConversion; + CVoigt_t C_voigt{CVoigt_t::Zero(vsize(DimM), vsize(DimM))}; + for (int i{0}; i < DimM; ++i) { + for (int j{0}; j < DimM; ++j) { + for (int k{0}; k < DimM; ++k) { + for (int l{0}; l < DimM; ++l) { + C_voigt(VC_t::sym_mat(i, j), VC_t::sym_mat(k, l)) = + get(C_col, i, j, k, l); + } + } + } + } + return C_voigt; + } + /* ---------------------------------------------------------------------- */ template std::vector MaterialBase::list_fields() const { return this->internal_fields.list_fields(); } template class MaterialBase<2, 2>; template class MaterialBase<2, 3>; template class MaterialBase<3, 3>; } // namespace muSpectre diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh index 6cd0285..e269c69 100644 --- a/src/materials/material_base.hh +++ b/src/materials/material_base.hh @@ -1,236 +1,245 @@ /** * @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 Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_MATERIALS_MATERIAL_BASE_HH_ #define SRC_MATERIALS_MATERIAL_BASE_HH_ #include "common/common.hh" #include "common/field.hh" #include "common/field_helpers.hh" #include "common/field_collection.hh" +#include "common/voigt_conversion.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; using LocalStressField_t = TensorField; //! Full type for strain fields using StrainField_t = StressField_t; using LocalStrainField_t = LocalStressField_t; //! Full type for tangent stiffness fields fields using TangentField_t = TensorField; using LocalTangentField_t = TensorField; using MScalarField_t = ScalarField, Real>; using MScalarFieldMap_t = ScalarFieldMap, Real, true>; using MVectorField_t = TensorField, Real, firstOrder, DimM>; using ScalarField_t = ScalarField, Real>; using VectorField_t = TensorField, Real, firstOrder, DimM>; + using CVoigt_t = Eigen::Matrix; + using Ccoord = Ccoord_t; //!< cell coordinates type using Rcoord = Rcoord_t; //!< real coordinates type //! Default constructor MaterialBase() = delete; //! Construct by name explicit 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 type need to overload add_pixel */ virtual void add_pixel(const Ccoord & ccooord); virtual void add_pixel_split(const Ccoord & local_ccoord, Real ratio); virtual void add_pixel_laminate(const Ccoord & local_ccoord, Real ratio, Rcoord normal_vec); // this function is responsible for allocating fields in case cells are // split or laminate void allocate_optional_fields(SplitCell is_cell_split = SplitCell::no); //! 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, SplitCell is_cell_split = SplitCell::no) = 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, SplitCell is_cell_split = SplitCell::no); //! computes stress and tangent moduli virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form, SplitCell is_cell_split = SplitCell::no) = 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, SplitCell is_cell_split = SplitCell::no); /** - * This functio gives the corresponding ratio of the material + * returns a voigt notation of stiffness matrix converted form column major + */ + + static auto make_C_voigt_from_col(T4Mat C_col) -> CVoigt_t; + + /** + * This functio gives the corresponding ratio of the material * assigned to the material */ Real get_assigned_ratio(Ccoord pixel); // this function return the ratio of which the // input pixel is consisted of this material void get_assigned_ratios(std::vector & pixel_assigned_ratios, Ccoord subdomain_resolutions, Ccoord subdomain_locations); //! 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(); } //! type to return real-valued fields in using EigenMap = Eigen::Map>; /** * return an internal field identified by its name as an Eigen Array */ EigenMap get_real_field(std::string field_name); /** * list the names of all internal fields */ std::vector list_fields() const; //! gives access to internal fields inline MFieldCollection_t & get_collection() { return this->internal_fields; } protected: const std::string name; //!< material's name (for output and debugging) MFieldCollection_t internal_fields{}; //!< storage for internal variables /** * //!< field holding the assigning ratio of the material * if the cell is split */ optional> assigned_ratio{}; /** *!< fields holding the: 1.normal vetor of the interface in laminate materials 2.local stress of the pixels 3 local stiffnesstangent materices of the pixels * if the cell is laminate */ optional> assigned_normal_vec{}; optional> local_stress_field{}; optional> local_tangent_field{}; private: }; } // namespace muSpectre #endif // SRC_MATERIALS_MATERIAL_BASE_HH_ diff --git a/src/materials/material_linear_elastic1.hh b/src/materials/material_linear_elastic1.hh index 4f6d234..cdc2505 100644 --- a/src/materials/material_linear_elastic1.hh +++ b/src/materials/material_linear_elastic1.hh @@ -1,188 +1,190 @@ /** * @file material_linear_elastic1.hh * * @author Till Junge * * @date 13 Nov 2017 * * @brief Implementation for linear elastic reference material like in de Geus * 2017. This follows the simplest and likely not most efficient * implementation (with exception of the Python law) * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC1_HH_ #define SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC1_HH_ #include "common/common.hh" #include "materials/stress_transformations_PK2.hh" #include "materials/material_muSpectre_base.hh" #include "materials/materials_toolbox.hh" namespace muSpectre { template class MaterialLinearElastic1; /** * traits for objective linear elasticity */ 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}; //! elasticity without internal variables using InternalVariables = std::tuple<>; }; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) /** * implements objective linear elasticity */ template class MaterialLinearElastic1 : public MaterialMuSpectre, DimS, DimM> { public: EIGEN_MAKE_ALIGNED_OPERATOR_NEW //! base class using Parent = MaterialMuSpectre; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename Parent::NeedTangent; //! global field collection 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; //! Hooke's law implementation using Hooke = typename MatTB::Hooke; //! Default constructor MaterialLinearElastic1() = delete; //! Copy constructor MaterialLinearElastic1(const MaterialLinearElastic1 & other) = delete; //! Construct by name, Young's modulus and Poisson's ratio MaterialLinearElastic1(std::string name, Real young, Real poisson); //! Move constructor MaterialLinearElastic1(MaterialLinearElastic1 && other) = delete; //! Destructor virtual ~MaterialLinearElastic1() = default; //! Copy assignment operator MaterialLinearElastic1 & operator=(const MaterialLinearElastic1 & other) = delete; //! Move assignment operator MaterialLinearElastic1 & operator=(MaterialLinearElastic1 && other) = delete; /** * evaluates second Piola-Kirchhoff stress given the Green-Lagrange * strain (or Cauchy stress if called with a small strain tensor) */ 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) */ template inline decltype(auto) evaluate_stress_tangent(s_t && E); /** * return the empty internals tuple */ InternalVariables & get_internals() { return this->internal_variables; } + auto get_C() { return this->C; } + protected: const Real young; //!< Young's modulus const Real poisson; //!< Poisson's ratio const Real lambda; //!< first Lamé constant const Real mu; //!< second Lamé constant (shear modulus) const Stiffness_t C; //!< stiffness tensor //! empty tuple InternalVariables internal_variables{}; private: }; /* ---------------------------------------------------------------------- */ template template auto MaterialLinearElastic1::evaluate_stress(s_t && E) -> decltype(auto) { return Hooke::evaluate_stress(this->lambda, this->mu, std::move(E)); } /* ---------------------------------------------------------------------- */ template template auto MaterialLinearElastic1::evaluate_stress_tangent(s_t && E) -> decltype(auto) { using Tangent_t = typename traits::TangentMap_t::reference; return Hooke::evaluate_stress( this->lambda, this->mu, Tangent_t(const_cast(this->C.data())), std::move(E)); } } // namespace muSpectre #endif // SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC1_HH_ diff --git a/src/materials/material_linear_elastic_generic1.cc b/src/materials/material_linear_elastic_generic1.cc index a82fdce..5fd00b6 100644 --- a/src/materials/material_linear_elastic_generic1.cc +++ b/src/materials/material_linear_elastic_generic1.cc @@ -1,71 +1,74 @@ /** * @file material_linear_elastic_generic1.cc * * @author Till Junge * * @date 21 Sep 2018 * * @brief implementation for MaterialLinearElasticGeneric * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "materials/material_linear_elastic_generic1.hh" #include "common/voigt_conversion.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ template MaterialLinearElasticGeneric1::MaterialLinearElasticGeneric1( const std::string & name, const CInput_t & C_voigt) - : Parent{name} { + : Parent{name} { using VC_t = VoigtConversion; constexpr Dim_t VSize{vsize(DimM)}; if (not(C_voigt.rows() == VSize) or not(C_voigt.cols() == VSize)) { std::stringstream err_str{}; err_str << "The stiffness tensor should be input as a " << VSize << " × " << VSize << " Matrix in Voigt notation. You supplied" << " a " << C_voigt.rows() << " × " << C_voigt.cols() << " matrix"; } for (int i{0}; i < DimM; ++i) { for (int j{0}; j < DimM; ++j) { for (int k{0}; k < DimM; ++k) { for (int l{0}; l < DimM; ++l) { get(this->C, i, j, k, l) = C_voigt(VC_t::sym_mat(i, j), VC_t::sym_mat(k, l)); } } } } } + /* ---------------------------------------------------------------------- */ + + template class MaterialLinearElasticGeneric1; template class MaterialLinearElasticGeneric1; template class MaterialLinearElasticGeneric1; } // namespace muSpectre diff --git a/src/materials/material_linear_elastic_generic1.hh b/src/materials/material_linear_elastic_generic1.hh index 37efcfb..dc851d5 100644 --- a/src/materials/material_linear_elastic_generic1.hh +++ b/src/materials/material_linear_elastic_generic1.hh @@ -1,186 +1,189 @@ /** * @file material_linear_elastic_generic1.hh * * @author Till Junge * * @date 21 Sep 2018 * * @brief Implementation fo a generic linear elastic material that * stores the full elastic stiffness tensor. Convenient but not the * most efficient * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC_GENERIC1_HH_ #define SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC_GENERIC1_HH_ #include "common/common.hh" #include "materials/stress_transformations_PK2.hh" #include "common/T4_map_proxy.hh" #include "materials/material_muSpectre_base.hh" #include "common/tensor_algebra.hh" namespace muSpectre { /** * forward declaration */ template class MaterialLinearElasticGeneric1; /** * traits for use by MaterialMuSpectre for crtp */ template struct MaterialMuSpectre_traits> { //! 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}; //! elasticity without internal variables using InternalVariables = std::tuple<>; }; /** * Linear elastic law defined by a full stiffness tensor. Very * generic, but not most efficient */ template class MaterialLinearElasticGeneric1 : public MaterialMuSpectre, DimS, DimM> { public: //! parent type using Parent = MaterialMuSpectre, DimS, DimM>; //! generic input tolerant to python input using CInput_t = Eigen::Ref, 0, Eigen::Stride>; + + using CVoigt_t = Eigen::Matrix; //! Default constructor MaterialLinearElasticGeneric1() = delete; /** * Constructor by name and stiffness tensor. * * @param name unique material name * @param C_voigt elastic tensor in Voigt notation */ MaterialLinearElasticGeneric1(const std::string & name, const CInput_t & C_voigt); //! Copy constructor MaterialLinearElasticGeneric1(const MaterialLinearElasticGeneric1 & other) = delete; //! Move constructor MaterialLinearElasticGeneric1(MaterialLinearElasticGeneric1 && other) = delete; //! Destructor virtual ~MaterialLinearElasticGeneric1() = default; //! Copy assignment operator MaterialLinearElasticGeneric1 & operator=(const MaterialLinearElasticGeneric1 & other) = delete; //! Move assignment operator MaterialLinearElasticGeneric1 & operator=(MaterialLinearElasticGeneric1 && other) = delete; //! see //! http://eigen.tuxfamily.org/dox/group__TopicStructHavingEigenMembers.html EIGEN_MAKE_ALIGNED_OPERATOR_NEW; /** * evaluates second Piola-Kirchhoff stress given the Green-Lagrange * strain (or Cauchy stress if called with a small strain tensor) */ template inline decltype(auto) evaluate_stress(const Eigen::MatrixBase & 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) */ template inline decltype(auto) evaluate_stress_tangent(const Eigen::MatrixBase & E); /** * return the empty internals tuple */ std::tuple<> & get_internals() { return this->internal_variables; } /** * return a reference to the stiffness tensor */ const T4Mat & get_C() const { return this->C; } + protected: T4Mat C{}; //! stiffness tensor //! empty tuple std::tuple<> internal_variables{}; }; /* ---------------------------------------------------------------------- */ template template auto MaterialLinearElasticGeneric1::evaluate_stress( const Eigen::MatrixBase & E) -> decltype(auto) { static_assert(Derived::ColsAtCompileTime == DimM, "wrong input size"); static_assert(Derived::RowsAtCompileTime == DimM, "wrong input size"); return Matrices::tensmult(this->C, E); } /* ---------------------------------------------------------------------- */ template template auto MaterialLinearElasticGeneric1::evaluate_stress_tangent( const Eigen::MatrixBase & E) -> decltype(auto) { using Stress_t = decltype(this->evaluate_stress(E)); using Stiffness_t = Eigen::Map>; using Ret_t = std::tuple; return Ret_t{this->evaluate_stress(E), Stiffness_t(this->C.data())}; } } // namespace muSpectre #endif // SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC_GENERIC1_HH_