diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp index ce18570..df742e1 100644 --- a/python/wrap/model.cpp +++ b/python/wrap/model.cpp @@ -1,259 +1,261 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas 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 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "model.hh" #include "adhesion_functional.hh" #include "functional.hh" #include "integral_operator.hh" #include "model_dumper.hh" #include "model_extensions.hh" #include "model_factory.hh" #include "numpy.hh" #include "residual.hh" #include "wrap.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ namespace wrap { using namespace py::literals; /// Wrap functional classes void wrapFunctionals(py::module& mod) { py::class_ func( mod, "Functional"); func.def(py::init<>()) .def("computeF", &functional::Functional::computeF) .def("computeGradF", &functional::Functional::computeGradF); py::class_ adh(mod, "AdhesionFunctional", func); adh.def_property("parameters", &functional::AdhesionFunctional::getParameters, &functional::AdhesionFunctional::setParameters) // legacy wrapper code .def("setParameters", &functional::AdhesionFunctional::setParameters); py::class_( mod, "ExponentialAdhesionFunctional", adh) .def(py::init&>(), "surface"_a); py::class_( mod, "MaugisAdhesionFunctional", adh) .def(py::init&>(), "surface"_a); py::class_( mod, "SquaredExponentialAdhesionFunctional", adh) .def(py::init&>(), "surface"_a); } template std::unique_ptr> instanciateFromNumpy(numpy& num) { std::unique_ptr> result = nullptr; switch (num.ndim()) { case 2: result = std::make_unique>>(num); return result; case 3: result = std::make_unique>>(num); return result; case 4: result = std::make_unique>>(num); return result; default: TAMAAS_EXCEPTION("instanciateFromNumpy expects the last dimension of numpy " "array to be the number of components"); } } /// Wrap IntegralOperator void wrapIntegralOperator(py::module& mod) { py::class_(mod, "IntegralOperator") .def("apply", [](IntegralOperator& op, numpy input, numpy output) { auto in = instanciateFromNumpy(input); auto out = instanciateFromNumpy(output); op.apply(*in, *out); }) .def("updateFromModel", &IntegralOperator::updateFromModel) .def("getModel", &IntegralOperator::getModel, py::return_value_policy::reference) .def("getKind", &IntegralOperator::getKind) .def("getType", &IntegralOperator::getType); } /// Wrap BEEngine classes void wrapBEEngine(py::module& mod) { py::class_(mod, "BEEngine") .def("solveNeumann", &BEEngine::solveNeumann) .def("solveDirichlet", &BEEngine::solveDirichlet) .def("getModel", &BEEngine::getModel, py::return_value_policy::reference); } /// Wrap Models void wrapModelClass(py::module& mod) { py::enum_(mod, "model_type") .value("basic_1d", model_type::basic_1d) .value("basic_2d", model_type::basic_2d) .value("surface_1d", model_type::surface_1d) .value("surface_2d", model_type::surface_2d) .value("volume_1d", model_type::volume_1d) .value("volume_2d", model_type::volume_2d); py::class_(mod, "Model") .def_property_readonly("type", &Model::getType) .def("setElasticity", &Model::setElasticity, "E"_a, "nu"_a) .def_property("E", &Model::getYoungModulus, &Model::setYoungModulus) .def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio) .def("getHertzModulus", &Model::getHertzModulus) .def("getYoungModulus", &Model::getYoungModulus) .def("getShearModulus", &Model::getShearModulus) .def("getPoissonRatio", &Model::getPoissonRatio) .def("getTraction", (GridBase & (Model::*)()) & Model::getTraction, py::return_value_policy::reference_internal) .def("getDisplacement", (GridBase & (Model::*)()) & Model::getDisplacement, py::return_value_policy::reference_internal) .def("getSystemSize", &Model::getSystemSize) .def("getDiscretization", &Model::getDiscretization) .def("getBoundarySystemSize", &Model::getBoundarySystemSize) .def("getBoundaryDiscretization", &Model::getBoundaryDiscretization) .def("solveNeumann", &Model::solveNeumann) .def("solveDirichlet", &Model::solveDirichlet) .def("dump", &Model::dump) .def("addDumper", &Model::addDumper, "dumper"_a) .def("setDumper", [](Model& m, const ModelDumper&) { TAMAAS_EXCEPTION("setDumper() is not a member of Model; use " "addDumper() instead"); }) .def("getBEEngine", &Model::getBEEngine, py::return_value_policy::reference_internal) .def("getIntegralOperator", &Model::getIntegralOperator, "operator_name"_a, py::return_value_policy::reference_internal) .def("registerField", [](Model& m, std::string name, numpy field) { auto f = instanciateFromNumpy(field); m.registerField(name, *f); }, "field_name"_a, "field"_a) .def("getField", &Model::getField, "field_name"_a, py::return_value_policy::reference_internal) .def("applyElasticity", [](Model& model, numpy stress, numpy strain) { auto out = instanciateFromNumpy(stress); auto in = instanciateFromNumpy(strain); model.applyElasticity(*out, *in); }) .def("__repr__", [](const Model& m) { std::stringstream ss; ss << m; return ss.str(); }) .def("__getitem__", [](const Model& m, std::string key) -> const GridBase& { try { return m.getField(key); } catch (std::out_of_range&) { throw py::key_error(); } }, py::return_value_policy::reference_internal); py::class_>( mod, "ModelDumper") .def(py::init<>()) .def("dump", &ModelDumper::dump, "model"_a) .def("__lshift__", [](ModelDumper& dumper, Model& model) { dumper << model; }); } /// Wrap factory for models void wrapModelFactory(py::module& mod) { py::class_(mod, "ModelFactory") .def_static("createModel", &ModelFactory::createModel, "model_type"_a, "system_size"_a, "discretization"_a) .def_static("createResidual", &ModelFactory::createResidual, - "model_type"_a, "hardening"_a, "sigma_y"_a); + "model_type"_a, "hardening"_a, "sigma_y"_a) + .def_static("registerVolumeOperators", + &ModelFactory::registerVolumeOperators, "model"_a); } /// Wrap residual class void wrapResidual(py::module& mod) { // TODO adapt to n-dim py::class_(mod, "Residual") .def(py::init()) .def("computeResidual", [](Residual& res, numpy& x) { auto in = instanciateFromNumpy(x); res.computeResidual(*in); }) .def("computeStress", [](Residual& res, numpy& x) { auto in = instanciateFromNumpy(x); res.computeStress(*in); }) .def("updateState", [](Residual& res, numpy& x) { auto in = instanciateFromNumpy(x); res.updateState(*in); }) .def("computeResidualDisplacement", [](Residual& res, numpy& x) { auto in = instanciateFromNumpy(x); res.computeResidualDisplacement(*in); }) .def("applyTangent", [](Residual& res, numpy& output, numpy& input, numpy& current_strain_inc) { auto out = instanciateFromNumpy(output); auto in = instanciateFromNumpy(input); auto inc = instanciateFromNumpy(current_strain_inc); res.applyTangent(*out, *in, *inc); }, "output"_a, "input"_a, "current_strain_increment"_a) .def("getVector", &Residual::getVector, py::return_value_policy::reference_internal) .def("getPlasticStrain", &Residual::getPlasticStrain, py::return_value_policy::reference_internal) .def("getStress", &Residual::getStress, py::return_value_policy::reference_internal); } void wrapModel(py::module& mod) { wrapBEEngine(mod); wrapModelClass(mod); wrapModelFactory(mod); wrapFunctionals(mod); wrapResidual(mod); wrapIntegralOperator(mod); } } // namespace wrap __END_TAMAAS__ diff --git a/src/model/elasto_plastic/residual.cpp b/src/model/elasto_plastic/residual.cpp index 71e6f14..a9b47b2 100644 --- a/src/model/elasto_plastic/residual.cpp +++ b/src/model/elasto_plastic/residual.cpp @@ -1,176 +1,172 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas 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 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "residual.hh" #include "grid_view.hh" +#include "model_factory.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ template void symmetrize(Grid& /*grad*/) { #if 0 Loop::loop( [](MatrixProxy gradient) { // Explicit loop here to avoid temporaries for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { gradient(i, j) = gradient(j, i) = 0.5 * (gradient(i, j) + gradient(j, i)); } } }, range>(grad)); #endif } /* -------------------------------------------------------------------------- */ template ResidualTemplate::ResidualTemplate(Model* model, Real sigma_0, Real h) : Residual(model), hardening(model, sigma_0, h) { - // Registering operators for residual calculation - model->registerIntegralOperator>("mindlin_gradient"); - model->registerIntegralOperator>("boussinesq_gradient"); - - // Registering operator for displacement calculation - model->registerIntegralOperator>("mindlin"); - model->registerIntegralOperator>("boussinesq"); + // Registering operators for residual and displacement calculation + ModelFactory::registerVolumeOperators(*model); // Initializing grids for (auto grid : {&strain, &stress, &residual, &tmp}) { grid->setNbComponents(trait::voigt); grid->resize(model->getDiscretization()); } // Registering fields model->registerField("stress", getStress()); model->registerField("plastic_strain", getPlasticStrain()); model->registerField("residual", getVector()); plastic_filter = [this](UInt layer) { #if 0 return true; #else return this->plastic_layers.find(layer) != this->plastic_layers.end(); #endif }; } /* -------------------------------------------------------------------------- */ template void ResidualTemplate::computeResidual(GridBase& strain_increment) { Grid& inc = dynamic_cast(strain_increment); hardening.template computePlasticIncrement(residual, strain, inc); this->updateFilter(residual); this->model->applyElasticity(residual, residual); integralOperator("mindlin_gradient") .applyIf(residual, residual, plastic_filter); integralOperator("boussinesq_gradient") .apply(this->model->getTraction(), tmp); residual -= strain_increment; residual += tmp; symmetrize(residual); } /* -------------------------------------------------------------------------- */ template void ResidualTemplate::applyTangent( GridBase& output, GridBase& input, GridBase& current_strain_increment) { Grid& inc = dynamic_cast(current_strain_increment); Grid& out = dynamic_cast(output); Grid& in = dynamic_cast(input); hardening.applyTangentIncrement(out, in, strain, inc); this->updateFilter(out); integralOperator("mindlin_gradient").applyIf(out, out, plastic_filter); out -= in; symmetrize(residual); } /* -------------------------------------------------------------------------- */ template void ResidualTemplate::computeStress(GridBase& strain_increment) { Grid& inc = dynamic_cast(strain_increment); hardening.template computeStress(stress, strain, inc); // integralOperator("boussinesq_gradient") // .apply(this->model->getTraction(), tmp); // this->model->applyElasticity(tmp, tmp); // stress += tmp; } /* -------------------------------------------------------------------------- */ template void ResidualTemplate::updateState( GridBase& converged_strain_increment) { Grid& inc = dynamic_cast(converged_strain_increment); hardening.template computeStress(stress, strain, inc); this->model->applyElasticity(residual, hardening.getPlasticStrain()); updateFilter(residual); strain += converged_strain_increment; // Computing displacements integralOperator("mindlin").applyIf(residual, this->model->getDisplacement(), plastic_filter); Grid disp_tmp(this->model->getDiscretization(), trait::components); integralOperator("boussinesq").apply(this->model->getTraction(), disp_tmp); this->model->getDisplacement() += disp_tmp; } /* -------------------------------------------------------------------------- */ template void ResidualTemplate::computeResidualDisplacement( GridBase& strain_increment) { Grid& inc = dynamic_cast(strain_increment); hardening.template computePlasticIncrement(tmp, strain, inc); updateFilter(tmp); this->model->applyElasticity(residual, tmp); integralOperator("mindlin").applyIf(residual, this->model->getDisplacement(), plastic_filter); } /* -------------------------------------------------------------------------- */ template void ResidualTemplate::updateFilter( Grid& plastic_strain_increment) { // plastic_layers.clear(); for (UInt i : Loop::range(plastic_strain_increment.sizes().front())) { auto slice = make_view(plastic_strain_increment, i); if (slice.dot(slice) / slice.dataSize() > 1e-14) plastic_layers.insert(i); } } /* -------------------------------------------------------------------------- */ template class ResidualTemplate; /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/src/model/model_factory.cpp b/src/model/model_factory.cpp index b43a8c8..bbbce11 100644 --- a/src/model/model_factory.cpp +++ b/src/model/model_factory.cpp @@ -1,49 +1,93 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas 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 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "model_factory.hh" #include "model_template.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ +#define MODEL_TYPE_CASE_MACRO(truc, arguments, type) \ + case type: { \ + result = std::make_unique>(std::forward(args)...); \ + break; \ + } + +#define IMPLEMENTED_MODEL_TYPES TAMAAS_MODEL_TYPES + +template class Concrete, + class... Args> +std::unique_ptr createFromModelType(model_type type, Args&&... args) { + std::unique_ptr result = nullptr; + + switch (type) { + BOOST_PP_SEQ_FOR_EACH(MODEL_TYPE_CASE_MACRO, ~, IMPLEMENTED_MODEL_TYPES); + default: + TAMAAS_EXCEPTION("Model type not implemented"); + } + + return result; +} + +#undef MODEL_TYPE_CASE_MACRO +#undef IMPLEMENTED_MODEL_TYPES + std::unique_ptr ModelFactory::createModel(model_type type, const std::vector& system_size, const std::vector& discretization) { - return createFromModelType(type, system_size, discretization); + return createFromModelType(type, system_size, + discretization); } std::unique_ptr ModelFactory::createResidual(Model* model, Real hardening, Real sigma_y) { if (model->getType() != model_type::volume_2d) TAMAAS_EXCEPTION("Cannot instanciate model: " << model); return std::make_unique>( model, hardening, sigma_y); } +#define MODEL_TYPE_CASE_MACRO(truc, arguments, type) \ + case type: { \ + m.registerIntegralOperator>("mindlin_gradient"); \ + m.registerIntegralOperator>("boussinesq_gradient"); \ + m.registerIntegralOperator>("mindlin"); \ + m.registerIntegralOperator>("boussinesq"); \ + break; \ + } + +void ModelFactory::registerVolumeOperators(Model& m) { + switch (m.getType()) { + BOOST_PP_SEQ_FOR_EACH(MODEL_TYPE_CASE_MACRO, ~, (model_type::volume_2d)); + default: + TAMAAS_EXCEPTION("Registering volume operators not supported on " << m); + } +} + +#undef MODEL_TYPE_CASE_MACRO __END_TAMAAS__ diff --git a/src/model/model_factory.hh b/src/model/model_factory.hh index f4bcbb3..4cae934 100644 --- a/src/model/model_factory.hh +++ b/src/model/model_factory.hh @@ -1,81 +1,58 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas 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 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __MODEL_FACTORY_HH__ #define __MODEL_FACTORY_HH__ /* -------------------------------------------------------------------------- */ #include "be_engine.hh" #include "model.hh" #include "model_type.hh" #include "residual.hh" #include #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /** * @brief Factory class for model * TODO put allocate grid in there (maybe) */ class ModelFactory { public: /// Create new model (TODO make std::unique_ptr work in swig) static std::unique_ptr createModel(model_type type, const std::vector& system_size, const std::vector& discretization); static std::unique_ptr createResidual(Model* model, Real hardening, Real sigma_y); + static void registerVolumeOperators(Model& m); }; -#define MODEL_TYPE_CASE_MACRO(truc, arguments, type) \ - case type: { \ - result = std::make_unique>(std::forward(args)...); \ - break; \ - } - -#define IMPLEMENTED_MODEL_TYPES TAMAAS_MODEL_TYPES - -template class Concrete, - class... Args> -std::unique_ptr createFromModelType(model_type type, Args&&... args) { - std::unique_ptr result = nullptr; - - switch (type) { - BOOST_PP_SEQ_FOR_EACH(MODEL_TYPE_CASE_MACRO, ~, IMPLEMENTED_MODEL_TYPES); - default: - TAMAAS_EXCEPTION("Model type not implemented"); - } - - return result; -} - -#undef MODEL_TYPE_CASE_MACRO -#undef IMPLEMENTED_MODEL_TYPES __END_TAMAAS__ #endif // __MODEL_FACTORY_HH__