diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp index da9df91..7131198 100644 --- a/python/wrap/model.cpp +++ b/python/wrap/model.cpp @@ -1,542 +1,520 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2022 Lucas Frérot * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program 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 Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "model.hh" #include "adhesion_functional.hh" #include "elastic_functional.hh" #include "field_container.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 #include #include /* -------------------------------------------------------------------------- */ namespace tamaas { namespace wrap { using namespace py::literals; struct model_operator_accessor { Model& m; decltype(auto) get(const std::string& name) { return m.getIntegralOperator(name); } }; 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"); } } void wrapFieldContainer(py::module& mod) { py::class_(mod, "FieldContainer") .def(py::init<>()) .def("__getitem__", &FieldContainer::at, py::return_value_policy::reference_internal) .def( "__setitem__", [](FieldContainer& f, std::string key, numpy grid) { f[key] = instanciateFromNumpy(grid); }, py::keep_alive<1, 3>()) .def( - "__setitem__", - [](FieldContainer& f, std::string key, numpy grid) { - f[key] = instanciateFromNumpy(grid); + "__contains__", + [](const FieldContainer& m, std::string key) { + const auto fields = m.fields(); + return std::find(fields.begin(), fields.end(), key) != fields.end(); + }, + "Test field existence") + .def( + "__iter__", + [](const FieldContainer& m) { + const auto& fields = m.fields_map(); + return py::make_key_iterator(fields.cbegin(), fields.cend()); + }, + py::keep_alive<0, 1>(), "Iterator on fields") + /// Deprecated interface + .def( + "registerField", + [](Model& m, std::string name, numpy field) { + TAMAAS_DEPRECATE("registerField()", "the [] operator"); + m[name] = instanciateFromNumpy(field); + }, + "field_name"_a, "field"_a, py::keep_alive<1, 3>()) + .def( + "getField", + [](const Model& m, std::string name) -> decltype(m.at(name)) { + TAMAAS_DEPRECATE("getField()", "the [] operator"); + return m.at(name); }, - py::keep_alive<1, 3>()); + "field_name"_a, py::return_value_policy::reference_internal) + .def( + "getFields", + [](const Model& m) { + TAMAAS_DEPRECATE("getFields()", "list(model)"); + return m.fields(); + }, + "Return fields list"); } /// Wrap functional classes void wrapFunctionals(py::module& mod) { py::class_, functional::wrap::PyFunctional> func(mod, "Functional"); func.def(py::init<>()) .def("computeF", &functional::Functional::computeF, "Compute functional value") .def("computeGradF", &functional::Functional::computeGradF, "Compute functional gradient"); py::class_( mod, "ElasticFunctionalPressure", func) .def(py::init&>()); py::class_(mod, "ElasticFunctionalGap", func) .def(py::init&>()); py::class_ adh(mod, "AdhesionFunctional", func); adh.def_property("parameters", &functional::AdhesionFunctional::getParameters, &functional::AdhesionFunctional::setParameters, "Parameters dictionary") .def("setParameters", [](functional::AdhesionFunctional& f, const std::map& m) { TAMAAS_DEPRECATE("setParameters()", "the parameters property"); f.setParameters(m); }); py::class_( mod, "ExponentialAdhesionFunctional", adh, "Potential of the form F = -γ·exp(-g/ρ)") .def(py::init&>(), "surface"_a); py::class_( mod, "MaugisAdhesionFunctional", adh, "Cohesive zone potential F = H(g - ρ)·γ/ρ") .def(py::init&>(), "surface"_a); py::class_( mod, "SquaredExponentialAdhesionFunctional", adh, "Potential of the form F = -γ·exp(-0.5·(g/ρ)²)") .def(py::init&>(), "surface"_a); } /// Wrap IntegralOperator void wrapIntegralOperator(py::module& mod) { py::class_(mod, "IntegralOperator") .def("apply", [](IntegralOperator& op, numpy input, numpy output) { TAMAAS_DEPRECATE("apply()", "the () operator"); auto in = instanciateFromNumpy(input); auto out = instanciateFromNumpy(output); op.apply(*in, *out); }) .def(TAMAAS_DEPRECATE_ACCESSOR(getModel, IntegralOperator, "model"), py::return_value_policy::reference) .def(TAMAAS_DEPRECATE_ACCESSOR(getKind, IntegralOperator, "kind")) .def(TAMAAS_DEPRECATE_ACCESSOR(getType, IntegralOperator, "type")) .def( "__call__", [](IntegralOperator& op, numpy input, numpy output) { auto in = instanciateFromNumpy(input); auto out = instanciateFromNumpy(output); op.apply(*in, *out); }, "Apply the integral operator") .def("updateFromModel", &IntegralOperator::updateFromModel, "Resets internal persistent variables from the model") .def_property_readonly("kind", &IntegralOperator::getKind) .def_property_readonly("model", &IntegralOperator::getModel) .def_property_readonly("type", &IntegralOperator::getType) .def_property_readonly("shape", &IntegralOperator::matvecShape) .def( "matvec", [](const IntegralOperator& op, numpy X) -> GridBase { GridBaseNumpy x(X); auto y = op.matvec(x); return y; }, py::return_value_policy::move); py::enum_(mod, "integration_method", "Integration method used for the computation " "of volumetric Fourier operators") .value("linear", integration_method::linear, "No approximation error, O(N₁·N₂·N₃) time complexity, may cause " "float overflow/underflow") .value("cutoff", integration_method::cutoff, "Approximation, O(sqrt(N₁²+N₂²)·N₃²) time complexity, no " "overflow/underflow risk"); } /// Wrap BEEngine classes void wrapBEEngine(py::module& mod) { py::class_(mod, "BEEngine") .def("solveNeumann", &BEEngine::solveNeumann) .def("solveDirichlet", &BEEngine::solveDirichlet) .def("registerNeumann", &BEEngine::registerNeumann) .def("registerDirichlet", &BEEngine::registerDirichlet) .def(TAMAAS_DEPRECATE_ACCESSOR(getModel, BEEngine, "model"), py::return_value_policy::reference) .def_property_readonly("model", &BEEngine::getModel); } template void wrapModelTypeTrait(py::module& mod) { using trait = model_type_traits; py::class_(mod, trait::repr) .def_property_readonly_static( "dimension", [](const py::object&) { return trait::dimension; }, "Dimension of computational domain") .def_property_readonly_static( "components", [](const py::object&) { return trait::components; }, "Number of components of vector fields") .def_property_readonly_static( "boundary_dimension", [](const py::object&) { return trait::boundary_dimension; }, "Dimension of boundary of computational domain") .def_property_readonly_static( "voigt", [](const py::object&) { return trait::voigt; }, "Number of components of symmetrical tensor fields") .def_property_readonly_static( "indices", [](const py::object&) { return trait::indices; }); } /// Wrap Models void wrapModelClass(py::module& mod) { py::enum_(mod, "model_type") .value("basic_1d", model_type::basic_1d, "Normal contact with 1D interface") .value("basic_2d", model_type::basic_2d, "Normal contact with 2D interface") .value("surface_1d", model_type::surface_1d, "Normal & tangential contact with 1D interface") .value("surface_2d", model_type::surface_2d, "Normal & tangential contact with 2D interface") .value("volume_1d", model_type::volume_1d, "Contact with volumetric representation and 1D interface") .value("volume_2d", model_type::volume_2d, "Contact with volumetric representation and 2D interface"); auto trait_mod = mod.def_submodule("_type_traits"); wrapModelTypeTrait(trait_mod); wrapModelTypeTrait(trait_mod); wrapModelTypeTrait(trait_mod); wrapModelTypeTrait(trait_mod); wrapModelTypeTrait(trait_mod); wrapModelTypeTrait(trait_mod); py::class_(mod, "_model_operator_acessor") .def(py::init()) .def( "__getitem__", [](model_operator_accessor& acc, std::string name) { try { return acc.get(name); } catch (std::out_of_range&) { throw py::key_error(name); } }, py::return_value_policy::reference_internal) .def("__contains__", [](model_operator_accessor& acc, std::string key) { const auto ops = acc.m.getIntegralOperators(); return std::find(ops.begin(), ops.end(), key) != ops.end(); }) .def( "__iter__", [](const model_operator_accessor& acc) { const auto& ops = acc.m.getIntegralOperatorsMap(); return py::make_key_iterator(ops.cbegin(), ops.cend()); }, py::keep_alive<0, 1>()); - py::class_(mod, "Model") + py::class_(mod, "Model") .def(py::init(py::overload_cast&, const std::vector&>( &ModelFactory::createModel))) .def_property_readonly("type", &Model::getType) .def_property("E", &Model::getYoungModulus, &Model::setYoungModulus, "Young's modulus") .def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio, "Poisson's ratio") .def_property_readonly("mu", &Model::getShearModulus, "Shear modulus") .def_property_readonly("E_star", &Model::getHertzModulus, "Contact (Hertz) modulus") .def_property_readonly("be_engine", &Model::getBEEngine, "Boundary element engine") .def( "setElasticity", [](Model& m, Real E, Real nu) { TAMAAS_DEPRECATE("setElasticity()", "the E and nu properties"); m.setElasticity(E, nu); }, "E"_a, "nu"_a) .def(TAMAAS_DEPRECATE_ACCESSOR(getHertzModulus, Model, "E_star")) .def(TAMAAS_DEPRECATE_ACCESSOR(getYoungModulus, Model, "E")) .def(TAMAAS_DEPRECATE_ACCESSOR(getShearModulus, Model, "mu")) .def(TAMAAS_DEPRECATE_ACCESSOR(getPoissonRatio, Model, "nu")) .def(TAMAAS_DEPRECATE_ACCESSOR(getTraction, Model, "traction"), py::return_value_policy::reference_internal) .def(TAMAAS_DEPRECATE_ACCESSOR(getDisplacement, Model, "displacement"), py::return_value_policy::reference_internal) .def(TAMAAS_DEPRECATE_ACCESSOR(getSystemSize, Model, "system_size")) .def(TAMAAS_DEPRECATE_ACCESSOR(getDiscretization, Model, "shape")) .def(TAMAAS_DEPRECATE_ACCESSOR(getBoundarySystemSize, Model, "boundary_system_size")) .def(TAMAAS_DEPRECATE_ACCESSOR(getBoundaryDiscretization, Model, "boundary_shape")) .def("solveNeumann", &Model::solveNeumann, "Solve surface tractions -> displacements") .def("solveDirichlet", &Model::solveDirichlet, "Solve surface displacemnts -> tractions") .def("dump", &Model::dump, "Write model data to registered dumpers") .def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>(), "Register a dumper") .def( "getBEEngine", [](Model& m) -> decltype(m.getBEEngine()) { TAMAAS_DEPRECATE("getBEEngine()", "the be_engine property"); return m.getBEEngine(); }, py::return_value_policy::reference_internal) .def( "getIntegralOperator", [](const Model& m, std::string name) { TAMAAS_DEPRECATE("getIntegralOperator()", "the operators property"); return m.getIntegralOperator(name); }, "operator_name"_a, py::return_value_policy::reference_internal) - .def( - "registerField", - [](Model& m, std::string name, numpy field) { - TAMAAS_DEPRECATE("registerField()", "the [] operator"); - m[name] = instanciateFromNumpy(field); - }, - "field_name"_a, "field"_a, py::keep_alive<1, 3>()) - .def( - "getField", - [](const Model& m, std::string name) -> decltype(m.at(name)) { - TAMAAS_DEPRECATE("getField()", "the [] operator"); - return m.at(name); - }, - "field_name"_a, py::return_value_policy::reference_internal) - .def( - "getFields", - [](const Model& m) { - TAMAAS_DEPRECATE("getFields()", "list(model)"); - return m.fields(); - }, - "Return fields list") - .def( "applyElasticity", [](Model& model, numpy stress, numpy strain) { auto out = instanciateFromNumpy(stress); auto in = instanciateFromNumpy(strain); model.applyElasticity(*out, *in); }, "Apply Hooke's law") // Python magic functions .def("__repr__", [](const Model& m) { std::stringstream ss; ss << m; return ss.str(); }) - .def( - "__getitem__", - [](const Model& m, std::string key) -> decltype(m.at(key)) { - try { - return m.at(key); - } catch (std::out_of_range&) { - throw py::key_error(key); - } - }, - py::return_value_policy::reference_internal, "Get field") - .def( - "__setitem__", - [](Model& m, std::string name, numpy field) { - m[name] = instanciateFromNumpy(field); - }, - py::keep_alive<1, 3>(), "Register new field") - .def( - "__contains__", - [](const Model& m, std::string key) { - const auto fields = m.fields(); - return std::find(fields.begin(), fields.end(), key) != fields.end(); - }, - py::keep_alive<0, 1>(), "Test field existence") - .def( - "__iter__", - [](const Model& m) { - const auto& fields = m.fields_map(); - return py::make_key_iterator(fields.cbegin(), fields.cend()); - }, - py::keep_alive<0, 1>(), "Iterator on fields") .def( "__copy__", [](const Model&) { throw std::runtime_error("__copy__ not implemented"); }, "Shallow copy of model. Not implemented.") .def( "__deepcopy__", [](const Model& m, py::dict) { return ModelFactory::createModel(m); }, "Deep copy of model.") .def_property_readonly("boundary_fields", &Model::getBoundaryFields) .def_property_readonly( "operators", [](Model& m) { return model_operator_accessor{m}; }, "Returns a dict-like object allowing access to the model's " "integral " "operators") // More python-like access to model properties .def_property_readonly("shape", &Model::getDiscretization, "Discretization (local in MPI environment)") .def_property_readonly("global_shape", &Model::getGlobalDiscretization, "Global discretization (in MPI environement)") .def_property_readonly("boundary_shape", &Model::getBoundaryDiscretization, "Number of points on boundary") .def_property_readonly("system_size", &Model::getSystemSize, "Size of physical domain") .def_property_readonly("boundary_system_size", &Model::getBoundarySystemSize, "Physical size of surface") .def_property_readonly("traction", (const GridBase& (Model::*)() const) & Model::getTraction, "Surface traction field") .def_property_readonly("displacement", (const GridBase& (Model::*)() const) & Model::getDisplacement, "Displacement field"); py::class_>( mod, "ModelDumper") .def(py::init<>()) .def("dump", &ModelDumper::dump, "model"_a, "Dump model") .def( "__lshift__", [](ModelDumper& dumper, Model& model) { dumper << model; }, "Dump model"); } /// Wrap factory for models void wrapModelFactory(py::module& mod) { py::class_(mod, "ModelFactory") .def_static( "createModel", py::overload_cast&, const std::vector&>( &ModelFactory::createModel), "model_type"_a, "system_size"_a, "global_discretization"_a, R"-(Create a new model of a given type, physical size and *global* discretization. :param model_type: the type of desired model :param system_size: the physical size of the domain in each direction :param global_discretization: number of points in each direction)-") .def_static("createModel", py::overload_cast(&ModelFactory::createModel), "model"_a, "Create a deep copy of a model.") .def_static("createResidual", &ModelFactory::createResidual, "model"_a, "sigma_y"_a, "hardening"_a = 0., R"-(Create an isotropic linear hardening residual. :param model: the model on which to define the residual :param sigma_y: the (von Mises) yield stress :param hardening: the hardening modulus)-") .def_static("registerVolumeOperators", &ModelFactory::registerVolumeOperators, "model"_a, "Register Boussinesq and Mindlin operators to model.") .def_static("setIntegrationMethod", &ModelFactory::setIntegrationMethod, "operator"_a, "method"_a, "cutoff"_a, "Set the integration method (linear or cutoff) for a volume " "integral operator"); } /// 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) .def("setIntegrationMethod", &Residual::setIntegrationMethod, "method"_a, "cutoff"_a = 1e-12) .def_property("yield_stress", &Residual::getYieldStress, &Residual::setYieldStress) .def_property("hardening_modulus", &Residual::getHardeningModulus, &Residual::setHardeningModulus) .def_property_readonly("model", &Residual::getModel); } void wrapModel(py::module& mod) { wrapFieldContainer(mod); wrapBEEngine(mod); wrapModelClass(mod); wrapModelFactory(mod); wrapFunctionals(mod); wrapResidual(mod); wrapIntegralOperator(mod); } } // namespace wrap } // namespace tamaas diff --git a/src/core/field_container.hh b/src/core/field_container.hh index ceb7c58..a477c34 100644 --- a/src/core/field_container.hh +++ b/src/core/field_container.hh @@ -1,76 +1,85 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2022 Lucas Frérot * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program 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 Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef FIELD_CONTAINER_HH #define FIELD_CONTAINER_HH /* -------------------------------------------------------------------------- */ #include "grid.hh" #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ -class FieldContainer { -public: +struct FieldContainer { + // Useful type names using Key = std::string; template using GridBasePtr = std::shared_ptr>; // TODO C++17: switch to std::variant using Value = boost::variant, GridBasePtr, GridBasePtr, GridBasePtr, GridBasePtr>; using FieldsMap = std::unordered_map; -public: + /// Destructor + virtual ~FieldContainer() = default; + + /// Access field pointer in const context const Value& at(const Key& name) const { return fields_.at(name); } + + /// Access/insert new pointer Value& operator[](const Key& name) { return fields_[name]; } + /// Access field of given type (non-const) template auto& field(const Key& name) { return *boost::get>(at(name)); } + /// Access field of given type (const) template const auto& field(const Key& name) const { return *boost::get>(at(name)); } + /// Return field keys decltype(auto) fields() const { std::vector keys; std::transform(fields_.cbegin(), fields_.cend(), std::back_inserter(keys), [](auto&& pair) { return pair.first; }); return keys; } + /// Return pointer map to fields const auto& fields_map() const { return fields_; } -protected: +private: FieldsMap fields_; }; } // namespace tamaas #endif diff --git a/src/model/model.cpp b/src/model/model.cpp index 84f2035..0a6a261 100644 --- a/src/model/model.cpp +++ b/src/model/model.cpp @@ -1,166 +1,166 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2022 Lucas Frérot * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program 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 Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "model.hh" #include "be_engine.hh" #include "logger.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ void Model::setElasticity(Real E_, Real nu_) { setYoungModulus(E_); setPoissonRatio(nu_); updateOperators(); } /* -------------------------------------------------------------------------- */ void Model::applyElasticity(GridBase& stress, const GridBase& strain) const { operators.at("hooke")->apply(const_cast&>(strain), stress); } /* -------------------------------------------------------------------------- */ GridBase& Model::getTraction() { return this->field("traction"); } const GridBase& Model::getTraction() const { return this->field("traction"); } /* -------------------------------------------------------------------------- */ GridBase& Model::getDisplacement() { return this->field("displacement"); } const GridBase& Model::getDisplacement() const { return this->field("displacement"); } /* -------------------------------------------------------------------------- */ const std::vector& Model::getSystemSize() const { return system_size; } /* -------------------------------------------------------------------------- */ const std::vector& Model::getDiscretization() const { return discretization; } /* -------------------------------------------------------------------------- */ void Model::solveNeumann() { engine->registerNeumann(); engine->solveNeumann(getTraction(), getDisplacement()); } void Model::solveDirichlet() { engine->registerDirichlet(); engine->solveDirichlet(getDisplacement(), getTraction()); } /* -------------------------------------------------------------------------- */ IntegralOperator* Model::getIntegralOperator(const std::string& name) const { return operators.at(name).get(); } /* -------------------------------------------------------------------------- */ std::vector Model::getIntegralOperators() const { std::vector keys; keys.reserve(operators.size()); std::transform(operators.begin(), operators.end(), std::back_inserter(keys), [](auto&& pair) { return pair.first; }); return keys; } /* -------------------------------------------------------------------------- */ void Model::updateOperators() { for (auto& op : operators) op.second->updateFromModel(); } /* -------------------------------------------------------------------------- */ void Model::addDumper(std::shared_ptr dumper) { this->dumpers.push_back(std::move(dumper)); } void Model::dump() const { for (const auto& dumper : dumpers) if (dumper) dumper->dump(*this); } /* -------------------------------------------------------------------------- */ std::vector Model::getBoundaryFields() const { std::vector boundary_fields; - for (auto&& pair : this->fields_) { + for (auto&& pair : this->fields_map()) { const auto boundary = boost::apply_visitor( [&](auto&& grid_ptr) { return this->isBoundaryField(*grid_ptr); }, pair.second); if (boundary) boundary_fields.push_back(pair.first); } return boundary_fields; } /* -------------------------------------------------------------------------- */ std::ostream& operator<<(std::ostream& o, const Model& _this) { o << "Model<" << _this.getType() << "> (E = " << _this.getYoungModulus() << ", nu = " << _this.getPoissonRatio() << ")\n"; auto out_collec = [&o](auto&& collec) { std::for_each(collec.begin(), collec.end() - 1, [&o](const auto& x) { o << x << ", "; }); o << collec.back(); }; // Printing domain size o << " - domain = ["; out_collec(_this.getSystemSize()); o << "]\n"; // Printing discretization o << " - discretization = ["; out_collec(_this.getDiscretization()); o << "]\n"; if (mpi::size() > 1) { o << " - global discretization = ["; out_collec(_this.getGlobalDiscretization()); o << "]\n"; } // Print fields o << " - registered fields = ["; out_collec(_this.fields()); o << "]\n"; o << " - registered operators = ["; out_collec(_this.getIntegralOperators()); o << "]"; if (not _this.dumpers.empty()) o << "\n - " << _this.dumpers.size() << " registered dumpers"; return o; } /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/src/model/model.hh b/src/model/model.hh index ca796cc..011935b 100644 --- a/src/model/model.hh +++ b/src/model/model.hh @@ -1,201 +1,197 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2022 Lucas Frérot * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Affero General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program 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 Affero General Public License for more details. * * You should have received a copy of the GNU Affero General Public License * along with this program. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef MODEL_HH #define MODEL_HH /* -------------------------------------------------------------------------- */ #include "be_engine.hh" #include "field_container.hh" #include "grid_base.hh" #include "integral_operator.hh" #include "model_dumper.hh" #include "model_type.hh" #include "tamaas.hh" #include #include #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * @brief Model containing pressure and displacement * This class is a container for the model fields. It is supposed to be * dimension agnostic, hence the GridBase members. */ class Model : public FieldContainer { protected: /// Constructor Model(std::vector system_size, std::vector discretization) - : system_size(std::move(system_size)), + : FieldContainer(), system_size(std::move(system_size)), discretization(std::move(discretization)) {} -public: - /// Destructor - virtual ~Model() = default; - public: /// Set elasticity parameters void setElasticity(Real E, Real nu); /// Get Hertz contact modulus Real getHertzModulus() const { return E / (1 - nu * nu); } /// Get Young's modulus Real getYoungModulus() const { return E; } /// Get Poisson's ratio Real getPoissonRatio() const { return nu; } /// Get shear modulus Real getShearModulus() const { return E / (2 * (1 + nu)); } /// Set Young's modulus void setYoungModulus(Real E_) { if (E_ < 0) TAMAAS_EXCEPTION("Elastic modulus should be positive"); this->E = E_; updateOperators(); } /// Set Poisson's ratio void setPoissonRatio(Real nu_) { if (nu_ > 0.5 or nu_ <= -1) TAMAAS_EXCEPTION("Poisson ratio should be in ]-1, 0.5]"); this->nu = nu_; updateOperators(); } // Model info accessors public: /// Get model type virtual model_type getType() const = 0; /// Get system physical size const std::vector& getSystemSize() const; /// Get boundary system physical size virtual std::vector getBoundarySystemSize() const = 0; /// Get discretization const std::vector& getDiscretization() const; /// Get discretization of global MPI system virtual std::vector getGlobalDiscretization() const = 0; /// Get boundary discretization virtual std::vector getBoundaryDiscretization() const = 0; /// Get boundary element engine BEEngine& getBEEngine() { TAMAAS_ASSERT(engine, "BEEngine was not initialized"); return *engine; } // Exposing some common operators public: /// Apply Hooke's law void applyElasticity(GridBase& stress, const GridBase& strain) const; /// Solve Neumann problem using default neumann operator void solveNeumann(); /// Solve Dirichlet problem using default dirichlet operator void solveDirichlet(); public: /// Register a new integral operator template IntegralOperator* registerIntegralOperator(const std::string& name); /// Get a registerd integral operator IntegralOperator* getIntegralOperator(const std::string& name) const; /// Get list of integral operators std::vector getIntegralOperators() const; /// Get operators mapcar const auto& getIntegralOperatorsMap() const { return operators; } /// Tell operators to update their cache void updateOperators(); public: /// Get pressure GridBase& getTraction(); /// Get pressure const GridBase& getTraction() const; /// Get displacement GridBase& getDisplacement(); /// Get displacement const GridBase& getDisplacement() const; /// Determine if a field is defined on boundary template bool isBoundaryField(const GridBase& field) const { return field.getDimension() == getBoundaryDiscretization().size(); } /// Return list of fields defined on boundary std::vector getBoundaryFields() const; public: /// Set the dumper object void addDumper(std::shared_ptr dumper); /// Dump the model void dump() const; friend std::ostream& operator<<(std::ostream& o, const Model& _this); protected: Real E = 1, nu = 0; std::vector system_size; std::vector discretization; std::unique_ptr engine = nullptr; std::unordered_map> operators; std::vector> dumpers; }; /* -------------------------------------------------------------------------- */ /* Template functions */ /* -------------------------------------------------------------------------- */ template IntegralOperator* Model::registerIntegralOperator(const std::string& name) { Logger().get(LogLevel::debug) << TAMAAS_DEBUG_MSG("registering operator " + name); operators[name] = std::make_unique(this); return operators[name].get(); } /* -------------------------------------------------------------------------- */ /* Output model to stream */ /* -------------------------------------------------------------------------- */ std::ostream& operator<<(std::ostream& o, const Model& _this); /* -------------------------------------------------------------------------- */ /* Simpler grid allocation */ /* -------------------------------------------------------------------------- */ template std::unique_ptr> allocateGrid(Model& model) { return allocateGrid( model.getType(), (boundary) ? model.getBoundaryDiscretization() : model.getDiscretization()); } template std::unique_ptr> allocateGrid(Model& model, UInt nc) { return allocateGrid(model.getType(), (boundary) ? model.getBoundaryDiscretization() : model.getDiscretization(), nc); } } // namespace tamaas #endif // MODEL_HH