diff --git a/python/tamaas/__init__.py.in b/python/tamaas/__init__.py.in index 646df6a..532da4c 100644 --- a/python/tamaas/__init__.py.in +++ b/python/tamaas/__init__.py.in @@ -1,83 +1,59 @@ # -*- mode:python; coding: utf-8 -*- # @file # @section LICENSE # # Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # 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 . """ Tamaas is a library dedicated to the fast treatment of rough contact problems. -See README.md and examples for guidance as how to use Tamaas. +See https://tamaas.readthedocs.io for documentation on how to use Tamaas. -See __authors__, __license__, __copyright__ for extra information about Tamaas. +See __author__, __license__, __copyright__ for extra information about Tamaas. """ __version__ = "@version@" __author__ = @authors@ # TODO Change copyright when is issue with unicode is found -__copyright__ = u"Copyright (©) 2016-20 EPFL " \ +__copyright__ = u"Copyright (©) 2016-2021 EPFL " \ + u"(École Polytechnique Fédérale de Lausanne), " \ + u"Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)" __license__ = "AGPLv3" __maintainer__ = "Lucas Frérot" __email__ = "@email@" try: from ._tamaas import model_type from ._tamaas import _type_traits as __tt - # Redefinition of model_type constants (for compatibility) - # mark as deprecated - model_type_basic_1d = model_type.basic_1d - model_type_basic_2d = model_type.basic_2d - model_type_surface_1d = model_type.surface_1d - model_type_surface_2d = model_type.surface_2d - model_type_volume_1d = model_type.volume_1d - model_type_volume_2d = model_type.volume_2d type_traits = { model_type.basic_1d: __tt.basic_1d, model_type.basic_2d: __tt.basic_2d, model_type.surface_1d: __tt.surface_1d, model_type.surface_2d: __tt.surface_2d, model_type.volume_1d: __tt.volume_1d, model_type.volume_2d: __tt.volume_2d, } del __tt from ._tamaas import * # noqa except ImportError as e: print("Error trying to import _tamaas:\n{}".format(e)) raise e - - -class ParamHelper: - """Legacy class to manage parameters/setters/getters""" - def __init__(self, obj): - self.obj = obj - - def set(self, params): - for key in params: - setter_name = 'set' + key - try: - accessor = getattr(self.obj, setter_name) - accessor(params[key]) - except Exception: - print("Setter '{}({})' does not exist for object {}" - .format(setter_name, type(params[key]), self.obj)) diff --git a/python/wrap.hh b/python/wrap.hh index 820fa33..d1d1d11 100644 --- a/python/wrap.hh +++ b/python/wrap.hh @@ -1,75 +1,91 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 __WRAP_HH__ #define __WRAP_HH__ /* -------------------------------------------------------------------------- */ -#include "tamaas.hh" -#include "numpy.hh" #include "cast.hh" +#include "numpy.hh" +#include "tamaas.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace tamaas { namespace wrap { namespace py = pybind11; /// For naming classes templated with dimension inline std::string makeDimensionName(const std::string& name, UInt dim) { std::stringstream str; str << name << dim << "D"; return str.str(); } -template class smart_pointer { +template +class smart_pointer { public: smart_pointer(T* ptr) : ptr(ptr) {} void assign(T val) { *ptr = val; } private: T* ptr; }; /* -------------------------------------------------------------------------- */ /* Prototypes for wrapping of tamaas components */ /* -------------------------------------------------------------------------- */ void wrapCore(py::module& mod); void wrapPercolation(py::module& mod); void wrapSurface(py::module& mod); void wrapModel(py::module& mod); void wrapSolvers(py::module& mod); void wrapTestFeatures(py::module& mod); void wrapCompute(py::module& mod); void wrapMPI(py::module& mod); #if defined(TAMAAS_LEGACY_BEM) void wrapBEM(py::module& mod); #endif +/* -------------------------------------------------------------------------- */ +/* Deprecation macro */ +/* -------------------------------------------------------------------------- */ +#define TAMAAS_DEPRECATE(olds, news) \ + do { \ + PyErr_WarnEx(PyExc_DeprecationWarning, \ + olds " is deprecated, use " news " instead.", 1); \ + } while (0) + +#define TAMAAS_DEPRECATE_ACCESSOR(acc, type, property) \ + #acc, [](const type& m) -> decltype(m.acc()) { \ + TAMAAS_DEPRECATE(#acc "()", "the " property " property"); \ + return m.acc(); \ + } + } // namespace wrap } // namespace tamaas #endif // __WRAP_HH__ diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp index 6d45377..6cfa48a 100644 --- a/python/wrap/model.cpp +++ b/python/wrap/model.cpp @@ -1,405 +1,431 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "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 /* -------------------------------------------------------------------------- */ namespace tamaas { namespace wrap { using namespace py::literals; struct model_operator_accessor { Model& m; decltype(auto) get(std::string name) { return m.getIntegralOperator(std::move(name)); } }; /// 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) .def("computeGradF", &functional::Functional::computeGradF); py::class_ adh(mod, "AdhesionFunctional", func); adh.def_property("parameters", &functional::AdhesionFunctional::getParameters, &functional::AdhesionFunctional::setParameters) - // TODO [legacy] remove wrapper code - .def("setParameters", &functional::AdhesionFunctional::setParameters); + .def("setParameters", + [](functional::AdhesionFunctional& f, + const std::map& m) -> void { + TAMAAS_DEPRECATE("setParameters()", "the parameters property"); + f.setParameters(m); + }); 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") - // TODO [legacy] remove .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("getModel", &IntegralOperator::getModel, + .def(TAMAAS_DEPRECATE_ACCESSOR(getModel, IntegralOperator, "model"), py::return_value_policy::reference) - .def("getKind", &IntegralOperator::getKind) - .def("getType", &IntegralOperator::getType) + .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); }) .def("updateFromModel", &IntegralOperator::updateFromModel) + .def_property_readonly("kind", &IntegralOperator::getKind) .def_property_readonly("model", &IntegralOperator::getModel) .def_property_readonly("type", &IntegralOperator::getType); py::enum_(mod, "integration_method") .value("linear", integration_method::linear) .value("cutoff", integration_method::cutoff); } /// 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) - // TODO [legacy] remove - .def("getModel", &BEEngine::getModel, py::return_value_policy::reference) + .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", [](py::object) { return trait::dimension; }) .def_property_readonly_static( "components", [](py::object) { return trait::components; }) .def_property_readonly_static( "boundary_dimension", [](py::object) { return trait::boundary_dimension; }) .def_property_readonly_static("voigt", [](py::object) { return trait::voigt; }) .def_property_readonly_static("indices", [](py::object) { return trait::indices; }); } /// 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); 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") .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) - // TODO [legacy] remove - .def("setElasticity", &Model::setElasticity, "E"_a, "nu"_a) - .def("getHertzModulus", &Model::getHertzModulus) - .def("getYoungModulus", &Model::getYoungModulus) - .def("getShearModulus", &Model::getShearModulus) - .def("getPoissonRatio", &Model::getPoissonRatio) - .def("getTraction", (GridBase & (Model::*)()) & Model::getTraction, + .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("getDisplacement", - (GridBase & (Model::*)()) & Model::getDisplacement, + .def(TAMAAS_DEPRECATE_ACCESSOR(getDisplacement, Model, "displacement"), py::return_value_policy::reference_internal) - .def("getSystemSize", &Model::getSystemSize) - .def("getDiscretization", &Model::getDiscretization) - .def("getBoundarySystemSize", &Model::getBoundarySystemSize) - .def("getBoundaryDiscretization", &Model::getBoundaryDiscretization) + .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) .def("solveDirichlet", &Model::solveDirichlet) .def("dump", &Model::dump) .def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>()) - .def("setDumper", - [](Model&, const ModelDumper&) { - TAMAAS_EXCEPTION("setDumper() is not a member of Model; use " - "addDumper() instead"); - }) - .def("getBEEngine", &Model::getBEEngine, + + .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", &Model::getIntegralOperator, + .def("getIntegralOperator", + [](const Model& m, std::string name) { + TAMAAS_DEPRECATE("getIntegralOperator()", + "the operators property"); + return m.getIntegralOperator(std::move(name)); + }, "operator_name"_a, py::return_value_policy::reference_internal) .def("registerField", [](Model& m, std::string name, numpy field) { + TAMAAS_DEPRECATE("registerField()", "the [] operator"); auto f = instanciateFromNumpy(field); m.registerField(name, std::move(f)); }, "field_name"_a, "field"_a, py::keep_alive<1, 3>()) .def("getField", - (GridBase & (Model::*)(const std::string&)) & Model::getField, + [](const Model& m, std::string name) -> decltype(m.getField(name)) { + TAMAAS_DEPRECATE("getField()", "the [] operator"); + return m.getField(std::move(name)); + }, "field_name"_a, py::return_value_policy::reference_internal) + .def("getFields", + [](const Model& m) { + TAMAAS_DEPRECATE("getFields()", "list(model)"); + return m.getFields(); + }, + "Return fields list") - .def("getFields", &Model::getFields, "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) -> const GridBase& { + [](const Model& m, std::string key) -> decltype(m[key]) { try { return m[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) { auto f = instanciateFromNumpy(field); m.registerField(name, std::move(f)); }, py::keep_alive<1, 3>(), "Register new field") .def("__contains__", [](const Model& m, std::string key) { const auto fields = m.getFields(); return std::find(fields.begin(), fields.end(), std::move(key)) != fields.end(); }, py::keep_alive<0, 1>(), "Test field existence") .def("__iter__", [](const Model& m) { const auto& fields = m.getFieldsMap(); return py::make_key_iterator(fields.cbegin(), fields.cend()); }, py::keep_alive<0, 1>(), "Iterator on fields") .def_property_readonly( "operators", [](Model& m) { return model_operator_accessor{m}; }, - "Returns a dict-like object allowing access to the model's integral " + "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", &ModelFactory::createModel, "model_type"_a, "system_size"_a, "global_discretization"_a, "Create a new model of a given type, physical size and " "*global* discretization") .def_static("createResidual", &ModelFactory::createResidual, "model"_a, "sigma_y"_a, "hardening"_a = 0., "Create an isotropic linear hardening residual") .def_static("registerVolumeOperators", &ModelFactory::registerVolumeOperators, "model"_a, "Register Boussinesq and Mindlin operators to model"); } /// 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) { wrapBEEngine(mod); wrapModelClass(mod); wrapModelFactory(mod); wrapFunctionals(mod); wrapResidual(mod); wrapIntegralOperator(mod); } } // namespace wrap } // namespace tamaas diff --git a/python/wrap/percolation.cpp b/python/wrap/percolation.cpp index c21ce0d..1fe3b26 100644 --- a/python/wrap/percolation.cpp +++ b/python/wrap/percolation.cpp @@ -1,96 +1,93 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "flood_fill.hh" #include "wrap.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace tamaas { namespace wrap { /* -------------------------------------------------------------------------- */ using namespace py::literals; /* -------------------------------------------------------------------------- */ template void wrapCluster(py::module& mod) { auto name = makeDimensionName("Cluster", dim); py::class_>(mod, name.c_str()) .def(py::init<>()) .def_property_readonly("area", &Cluster::getArea, "Area of cluster") .def_property_readonly("points", &Cluster::getPoints, "Get list of points of cluster") .def_property_readonly("perimeter", &Cluster::getPerimeter, "Get perimeter of cluster") .def_property_readonly("bounding_box", &Cluster::boundingBox, "Compute the bounding box of a cluster") - // TODO [legacy] remove - .def("getArea", &Cluster::getArea, "Area of cluster") - .def("getPoints", &Cluster::getPoints, - "Get list of points of cluster") - .def("getPerimeter", &Cluster::getPerimeter, - "Get perimeter of cluster") + .def(TAMAAS_DEPRECATE_ACCESSOR(getArea, Cluster, "area")) + .def(TAMAAS_DEPRECATE_ACCESSOR(getPoints, Cluster, "points")) + .def(TAMAAS_DEPRECATE_ACCESSOR(getPerimeter, Cluster, "perimeter")) .def("__str__", [](const Cluster& cluster) { std::stringstream sstr; sstr << '{'; auto&& points = cluster.getPoints(); auto print_point = [&sstr](const auto& point) { sstr << '('; for (UInt i = 0; i < point.size() - 1; ++i) sstr << point[i] << ", "; sstr << point.back() << ')'; }; auto point_it = points.begin(); for (UInt p = 0; p < points.size() - 1; ++p) { print_point(*(point_it++)); sstr << ", "; } print_point(points.back()); sstr << "}"; return sstr.str(); }); } void wrapPercolation(py::module& mod) { wrapCluster<1>(mod); wrapCluster<2>(mod); wrapCluster<3>(mod); py::class_(mod, "FloodFill") .def_static("getSegments", &FloodFill::getSegments, "Return a list of segments from boolean map", "contact"_a) .def_static("getClusters", &FloodFill::getClusters, "Return a list of clusters from boolean map", "contact"_a, "diagonal"_a) .def_static("getVolumes", &FloodFill::getVolumes, "Return a list of volume clusters", "map"_a, "diagonal"_a); } } // namespace wrap } // namespace tamaas diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp index e5c01bf..300b2b1 100644 --- a/python/wrap/solvers.cpp +++ b/python/wrap/solvers.cpp @@ -1,179 +1,188 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "beck_teboulle.hh" #include "condat.hh" #include "contact_solver.hh" #include "dfsane_solver.hh" #include "ep_solver.hh" #include "epic.hh" #include "kato.hh" #include "kato_saturated.hh" #include "polonsky_keer_rey.hh" #include "polonsky_keer_tan.hh" #include "wrap.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { namespace wrap { /* -------------------------------------------------------------------------- */ using namespace py::literals; class PyEPSolver : public EPSolver { public: using EPSolver::EPSolver; void solve() override { PYBIND11_OVERLOAD_PURE(void, EPSolver, solve); } void updateState() override { PYBIND11_OVERLOAD(void, EPSolver, updateState); } }; /* -------------------------------------------------------------------------- */ void wrapSolvers(py::module& mod) { py::class_(mod, "ContactSolver") - // TODO [legacy] remove - .def("setMaxIterations", &ContactSolver::setMaxIterations, "max_iter"_a) - .def("setDumpFrequency", &ContactSolver::setDumpFrequency, "dump_freq"_a) + .def("setMaxIterations", + [](ContactSolver& m, UInt iter) { + TAMAAS_DEPRECATE("setMaxIterations()", "the max_iter property"); + m.setMaxIterations(iter); + }, + "max_iter"_a) + .def("setDumpFrequency", + [](ContactSolver& m, UInt freq) { + TAMAAS_DEPRECATE("setDumpFrequency()", "the dump_freq property"); + m.setDumpFrequency(freq); + }, + "dump_freq"_a) .def_property("max_iter", &ContactSolver::getMaxIterations, &ContactSolver::setMaxIterations, "Maximum number of iterations") .def_property("dump_freq", &ContactSolver::getDumpFrequency, &ContactSolver::setDumpFrequency, "Frequency of displaying solver info") .def_property_readonly("model", &ContactSolver::getModel) .def_property_readonly("surface", &ContactSolver::getSurface) .def("addFunctionalTerm", &ContactSolver::addFunctionalTerm) .def("solve", py::overload_cast>(&ContactSolver::solve), "target_force"_a, py::call_guard()) .def("solve", py::overload_cast(&ContactSolver::solve), "target_normal_pressure"_a, py::call_guard()); py::class_ pkr(mod, "PolonskyKeerRey"); // Need to export enum values before defining PKR constructor py::enum_(pkr, "type") .value("gap", PolonskyKeerRey::gap) .value("pressure", PolonskyKeerRey::pressure) .export_values(); pkr.def(py::init&, Real, PolonskyKeerRey::type, PolonskyKeerRey::type>(), "model"_a, "surface"_a, "tolerance"_a, "primal_type"_a = PolonskyKeerRey::type::pressure, "constraint_type"_a = PolonskyKeerRey::type::pressure, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("computeError", &PolonskyKeerRey::computeError); py::class_(mod, "KatoSaturated") .def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "pmax"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def_property("pmax", &KatoSaturated::getPMax, &KatoSaturated::setPMax); py::class_ kato(mod, "Kato"); kato.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &Kato::solve, "p0"_a, "proj_iter"_a = 50) .def("solveRelaxed", &Kato::solveRelaxed, "g0"_a) .def("solveRegularized", &Kato::solveRegularized, "p0"_a, "r"_a = 0.01) .def("computeCost", &Kato::computeCost, "use_tresca"_a = false); py::class_ bt(mod, "BeckTeboulle"); bt.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &BeckTeboulle::solve, "p0"_a) .def("computeCost", &BeckTeboulle::computeCost); py::class_ cd(mod, "Condat"); cd.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &Condat::solve, "p0"_a, "grad_step"_a = 0.9) .def("computeCost", &Condat::computeCost); py::class_ pkt(mod, "PolonskyKeerTan"); pkt.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &PolonskyKeerTan::solve, "p0"_a) .def("solveTresca", &PolonskyKeerTan::solveTresca, "p0"_a) .def("computeCost", &PolonskyKeerTan::computeCost, "use_tresca"_a = false); py::class_(mod, "_tolerance_manager") .def(py::init()); py::class_(mod, "EPSolver") .def(py::init(), "residual"_a, py::keep_alive<1, 2>()) .def("solve", &EPSolver::solve) .def("getStrainIncrement", &EPSolver::getStrainIncrement, py::return_value_policy::reference_internal) .def("getResidual", &EPSolver::getResidual, py::return_value_policy::reference_internal) .def("updateState", &EPSolver::updateState) .def_property("tolerance", &EPSolver::getTolerance, &EPSolver::setTolerance) .def("setToleranceManager", &EPSolver::setToleranceManager) .def("beforeSolve", &EPSolver::beforeSolve); py::class_(mod, "_DFSANESolver") .def(py::init(), "residual"_a, py::keep_alive<1, 2>()) .def(py::init([](Residual& res, Model&) { return std::make_unique(res); }), "residual"_a, "model"_a, py::keep_alive<1, 2>(), "Deprecated constructor: only here for compatibility reasons"); py::class_(mod, "EPICSolver") .def(py::init(), "contact_solver"_a, "elasto_plastic_solver"_a, "tolerance"_a = 1e-10, "relaxation"_a = 0.3, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", [](EPICSolver& solver, Real pressure) { return solver.solve(std::vector{pressure}); }, "normal_pressure"_a, py::call_guard()) .def("acceleratedSolve", [](EPICSolver& solver, Real pressure) { return solver.acceleratedSolve(std::vector{pressure}); }, "normal_pressure"_a, py::call_guard()); } } // namespace wrap } // namespace tamaas diff --git a/python/wrap/surface.cpp b/python/wrap/surface.cpp index cbce617..f5bad72 100644 --- a/python/wrap/surface.cpp +++ b/python/wrap/surface.cpp @@ -1,217 +1,215 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "isopowerlaw.hh" #include "regularized_powerlaw.hh" #include "surface_generator_filter.hh" #include "surface_generator_filter_fft.hh" #include "surface_generator_random_phase.hh" #include "wrap.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ namespace wrap { using namespace py::literals; /* -------------------------------------------------------------------------- */ template class PyFilter : public Filter { public: using Filter::Filter; // Overriding pure virtual functions void computeFilter(GridHermitian& filter_coefficients) const override { PYBIND11_OVERLOAD_PURE(void, Filter, computeFilter, filter_coefficients); } }; template void wrapFilter(py::module& mod) { auto name = makeDimensionName("Filter", dim); py::class_, std::shared_ptr>, PyFilter>( mod, name.c_str()) .def(py::init<>()) .def("computeFilter", (void (Filter::*)(GridHermitian&) const) & Filter::computeFilter); } /* -------------------------------------------------------------------------- */ template void wrapIsopowerlaw(py::module& mod) { std::string name = makeDimensionName("Isopowerlaw", dim); py::class_, Filter, std::shared_ptr>>( mod, name.c_str()) .def(py::init<>()) .def_property("q0", &Isopowerlaw::getQ0, &Isopowerlaw::setQ0, "Long wavelength cutoff") .def_property("q1", &Isopowerlaw::getQ1, &Isopowerlaw::setQ1, "Rolloff wavelength") .def_property("q2", &Isopowerlaw::getQ2, &Isopowerlaw::setQ2, "Short wavelength cutoff") .def_property("hurst", &Isopowerlaw::getHurst, &Isopowerlaw::setHurst, "Hurst exponent") .def("rmsHeights", &Isopowerlaw::rmsHeights, "Theoretical RMS of heights") .def("moments", &Isopowerlaw::moments, "Theoretical first 3 moments of spectrum") .def("alpha", &Isopowerlaw::alpha, "Nayak's bandwidth parameter") .def("rmsSlopes", &Isopowerlaw::rmsSlopes, - "Theoretical RMS of slopes") - // TODO [legacy] remove wrapper code - .def("getQ0", &Isopowerlaw::getQ0, - py::return_value_policy::reference) - .def("getQ1", &Isopowerlaw::getQ1, - py::return_value_policy::reference) - .def("getQ2", &Isopowerlaw::getQ2, - py::return_value_policy::reference) - .def("setQ0", &Isopowerlaw::setQ0, - py::return_value_policy::reference) - .def("setQ1", &Isopowerlaw::setQ1, - py::return_value_policy::reference) - .def("setQ2", &Isopowerlaw::setQ2, - py::return_value_policy::reference) - .def("setHurst", &Isopowerlaw::setHurst) - .def("getHurst", &Isopowerlaw::getHurst, - py::return_value_policy::reference); + "Theoretical RMS of slopes"); name = makeDimensionName("RegularizedPowerlaw", dim); py::class_, Filter, std::shared_ptr>>(mod, name.c_str()) .def(py::init<>()) .def_property("q1", &RegularizedPowerlaw::getQ1, &RegularizedPowerlaw::setQ1, "Long wavelength cutoff") .def_property("q2", &RegularizedPowerlaw::getQ2, &RegularizedPowerlaw::setQ2, "Short wavelength cutoff") .def_property("hurst", &RegularizedPowerlaw::getHurst, &RegularizedPowerlaw::setHurst, "Hurst exponent"); } /* -------------------------------------------------------------------------- */ template void wrapSurfaceGenerators(py::module& mod) { std::string generator_name = makeDimensionName("SurfaceGenerator", dim); py::class_>(mod, generator_name.c_str()) .def("buildSurface", &SurfaceGenerator::buildSurface, py::return_value_policy::reference_internal) + .def("setSizes", - py::overload_cast>( - &SurfaceGenerator::setSizes), - "global_surface_shape"_a) + [](SurfaceGenerator& m, std::array s) { + TAMAAS_DEPRECATE("setSizes()", "the shape property"); + m.setSizes(std::move(s)); + }) + .def("setRandomSeed", + [](SurfaceGenerator& m, long s) { + TAMAAS_DEPRECATE("setRandomSeed()", "the random_seed property"); + m.setRandomSeed(s); + }) + .def_property("random_seed", &SurfaceGenerator::getRandomSeed, &SurfaceGenerator::setRandomSeed, "Random generator seed") .def_property("shape", &SurfaceGenerator::getSizes, py::overload_cast>( &SurfaceGenerator::setSizes), "Global shape of surfaces"); std::string filter_name = makeDimensionName("SurfaceGeneratorFilter", dim); py::class_, SurfaceGenerator>( mod, filter_name.c_str()) .def(py::init<>(), "Default constructor") .def(py::init>(), "Initialize with global surface shape") - // TODO [legacy] remove - .def("setFilter", &SurfaceGeneratorFilter::setFilter, + + .def("setFilter", + [](SurfaceGeneratorFilter& m, std::shared_ptr> s) { + TAMAAS_DEPRECATE("setFilter()", "the spectrum property"); + m.setSpectrum(std::move(s)); + }, + "Set PSD filter", "filter"_a) + .def("setSpectrum", + [](SurfaceGeneratorFilter& m, std::shared_ptr> s) { + TAMAAS_DEPRECATE("setSpectrum()", "the spectrum property"); + m.setSpectrum(std::move(s)); + }, "Set PSD filter", "filter"_a) - .def("setSpectrum", &SurfaceGeneratorFilter::setSpectrum, - "Set PSD spectrum", "filter"_a) - .def("setRandomSeed", &SurfaceGeneratorFilter::setRandomSeed, - "[[Deprecated: use property]]") .def_property("spectrum", &SurfaceGeneratorFilter::getSpectrum, &SurfaceGeneratorFilter::setSpectrum, "Power spectrum object"); std::string random_name = makeDimensionName("SurfaceGeneratorRandomPhase", dim); py::class_, SurfaceGeneratorFilter>( mod, random_name.c_str()) .def(py::init<>(), "Default constructor") .def(py::init>(), "Initialize with global surface shape"); } /* -------------------------------------------------------------------------- */ // TODO [legacy] remove wrap void wrapSurfaceGeneratorFilterFFT(py::module& mod) { py::class_(mod, "SurfaceGeneratorFilterFFT") .def(py::init<>()) .def("getQ0", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getQ0()); }) .def("getQ1", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getQ1()); }) .def("getQ2", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getQ2()); }) .def("getHurst", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getHurst()); }) .def("analyticalRMS", &SurfaceGeneratorFilterFFT::analyticalRMS) .def("buildSurface", &SurfaceGeneratorFilterFFT::buildSurface, py::return_value_policy::reference_internal) .def("getGridSize", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getGridSize()); }) .def("getRandomSeed", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getRandomSeed()); }) .def("getRMS", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getRMS()); }) .def("Init", &SurfaceGeneratorFilterFFT::Init); } /* -------------------------------------------------------------------------- */ void wrapSurface(py::module& mod) { wrapFilter<1>(mod); wrapFilter<2>(mod); wrapIsopowerlaw<1>(mod); wrapIsopowerlaw<2>(mod); wrapSurfaceGenerators<1>(mod); wrapSurfaceGenerators<2>(mod); // legacy wrap #if defined(TAMAAS_LEGACY_BEM) wrapSurfaceGeneratorFilterFFT(mod); #endif } } // namespace wrap /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/src/model/model.cpp b/src/model/model.cpp index eb5df5f..27d435e 100644 --- a/src/model/model.cpp +++ b/src/model/model.cpp @@ -1,184 +1,184 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 getField("traction"); } const GridBase& Model::getTraction() const { return (*this)["traction"]; } /* -------------------------------------------------------------------------- */ GridBase& Model::getDisplacement() { return getField("displacement"); } const GridBase& Model::getDisplacement() const { return (*this)["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(getTraction(), getDisplacement()); } /* -------------------------------------------------------------------------- */ -IntegralOperator* Model::getIntegralOperator(const std::string& name) { +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 (auto& dumper : dumpers) if (dumper) dumper->dump(*this); } /* -------------------------------------------------------------------------- */ void Model::registerField(const std::string& name, std::shared_ptr> field) { fields[name] = std::move(field); } const GridBase& Model::getField(const std::string& name) const try { return *fields.at(name); } catch (std::out_of_range& e) { Logger().get(LogLevel::warning) << "Field " << name << " not registered in model\n"; throw e; } GridBase& Model::getField(const std::string& name) try { return *fields.at(name); } catch (std::out_of_range& e) { Logger().get(LogLevel::warning) << "Field " << name << " not registered in model\n"; throw e; } std::vector Model::getFields() const { std::vector keys; keys.reserve(fields.size()); std::transform(fields.begin(), fields.end(), std::back_inserter(keys), [](auto&& pair) { return pair.first; }); return keys; } const GridBase& Model::operator[](const std::string& name) const { return getField(name); } GridBase& Model::operator[](const std::string& name) { return getField(name); } /* -------------------------------------------------------------------------- */ 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.getFields()); o << "]\n"; o << " - registered operators = ["; out_collec(_this.getIntegralOperators()); o << "]"; if (_this.dumpers.size()) o << "\n - " << _this.dumpers.size() << " registered dumpers"; return o; } /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/src/model/model.hh b/src/model/model.hh index b454ee6..687dcb7 100644 --- a/src/model/model.hh +++ b/src/model/model.hh @@ -1,210 +1,210 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "grid_base.hh" #include "integral_operator.hh" #include "model_dumper.hh" #include "model_type.hh" #include "tamaas.hh" #include #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 { protected: /// Constructor Model(std::vector system_size, std::vector discretization) : 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); + 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; /// Register a field void registerField(const std::string& name, std::shared_ptr> field); /// Get a field const GridBase& getField(const std::string& name) const; /// Get a non-const field GridBase& getField(const std::string& name); /// Get fields std::vector getFields() const; /// Get fields map const auto& getFieldsMap() const { return fields; } /// Get a field const GridBase& operator[](const std::string& name) const; /// Get a non-const field GridBase& operator[](const std::string& name); 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::unordered_map>> fields; 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 diff --git a/src/percolation/flood_fill.hh b/src/percolation/flood_fill.hh index e5ad746..3f3554e 100644 --- a/src/percolation/flood_fill.hh +++ b/src/percolation/flood_fill.hh @@ -1,77 +1,77 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 FLOOD_FILL_H #define FLOOD_FILL_H /* -------------------------------------------------------------------------- */ #include "grid.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ template class Cluster { using Point = std::array; public: /// Constructor Cluster(Point start, const Grid& map, Grid& visited, bool diagonal); /// Copy constructor Cluster(const Cluster& other); /// Default constructor Cluster() = default; /// Get area of cluster - UInt getArea() { return getPoints().size(); } + UInt getArea() const { return getPoints().size(); } /// Get perimeter of cluster - UInt getPerimeter() { return perimeter; } + UInt getPerimeter() const { return perimeter; } /// Get contact points const std::list& getPoints() const { return points; } /// Get bounding box std::pair, std::array> boundingBox() const; /// Assign next neighbors auto getNextNeighbors(const std::array& p); /// Assign diagonal neighbors auto getDiagonalNeighbors(const std::array& p); private: std::list points; UInt perimeter = 0; }; /* -------------------------------------------------------------------------- */ class FloodFill { public: static std::list> getSegments(const Grid& contact); /// Returns a list of clusters from a boolean contact map static std::list> getClusters(const Grid& contact, bool diagonal); static std::list> getVolumes(const Grid& map, bool diagonal); }; /* -------------------------------------------------------------------------- */ } // namespace tamaas /* -------------------------------------------------------------------------- */ #endif // FLOOD_FILL_H diff --git a/tests/test_boussinesq_surface.py b/tests/test_boussinesq_surface.py index 180076d..9cb9333 100644 --- a/tests/test_boussinesq_surface.py +++ b/tests/test_boussinesq_surface.py @@ -1,80 +1,80 @@ # -*- coding: utf-8 -*- # @file # @section LICENSE # # Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # 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 . import numpy as np import tamaas as tm from numpy.linalg import norm def test_boussinesq_surface(tamaas_fixture): # Definition of modeled domain model_type = tm.model_type.volume_2d discretization = [2, 128, 128] system_size = [1., 1., 1.] x, y = [np.linspace(0, size, n, endpoint=False, dtype=tm.dtype) for size, n in zip(system_size[1:], discretization[1:])] xx, yy = np.meshgrid(x, y, sparse=True) tractions = np.zeros(discretization[1:] + [3], dtype=tm.dtype) for i in range(3): omega = np.random.randint(1, 20, (2,)) * 2 * np.pi tractions[:, :, i] = np.cos(omega[0]*xx)*np.cos(omega[1]*yy) # Material contants E = 1. # Young's modulus nu = 0.3 # Poisson's ratio # Creation of model model_vol = tm.ModelFactory.createModel(model_type, system_size, discretization) model_surf = tm.ModelFactory.createModel(tm.model_type.surface_2d, [1, 1], tractions.shape[:-1]) model_vol.E = model_surf.E = E model_vol.nu = model_surf.nu = nu # Setup for integral operators tm.ModelFactory.createResidual(model_vol, 0, 0) # Pressure definition model_vol['traction'][...] = tractions model_surf['traction'][...] = tractions model_surf.solveNeumann() # Applying operator boussinesq = model_vol.operators["boussinesq"] boussinesq(model_vol['traction'], model_vol['displacement']) # # Dumper # dumper_helper = UVWDumper(residual, args.name) # model.addDumper(dumper_helper) # model.dump() # print("Done") for i in range(3): - error = norm(model_vol.getDisplacement()[0, :, :, i] - - model_surf.getDisplacement()[:, :, i]) \ + error = norm(model_vol.displacement[0, :, :, i] + - model_surf.displacement[:, :, i]) \ / norm(tractions[:, :, i]) assert error < 1e-16 diff --git a/tests/test_hertz_disp.py b/tests/test_hertz_disp.py index 19db312..3ff22ee 100644 --- a/tests/test_hertz_disp.py +++ b/tests/test_hertz_disp.py @@ -1,133 +1,133 @@ # -*- coding: utf-8 -*- # @file # @section LICENSE # # Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # 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 . from __future__ import print_function, division import tamaas as tm import numpy as np def constructHertzProfile(size, curvature): radius = 1. / curvature x = np.linspace(-0.5, 0.5, size, dtype=tm.dtype) y = np.linspace(-0.5, 0.5, size, dtype=tm.dtype) x, y = np.meshgrid(x, y) surface = np.sqrt(radius**2 - x**2 - y**2) surface -= surface.min() return surface.copy() def computeHertzDisplacement(e_star, contact_size, max_pressure, size): x = np.linspace(-0.5, 0.5, size) y = np.linspace(-0.5, 0.5, size) x, y = np.meshgrid(x, y) disp = np.pi * max_pressure / (4 * contact_size * e_star) \ * (2 * contact_size**2 - (x**2 + y**2)) return disp.copy() def test_hertz_disp(tamaas_fixture): grid_size = 512 curvature = 0.5 effective_modulus = 1. load = 0.12 surface = constructHertzProfile(grid_size, curvature) - model = tm.ModelFactory.createModel(tm.model_type_basic_2d, + model = tm.ModelFactory.createModel(tm.model_type.basic_2d, [1., 1.], [grid_size, grid_size]) model.E, model.nu = 1, 0 solver = tm.PolonskyKeerRey(model, surface, 1e-12, tm.PolonskyKeerRey.gap, tm.PolonskyKeerRey.gap) solver.solve(load - surface.mean()) tractions = model['traction'] true_displacements = model['displacement'] true_gap = true_displacements - surface pressure = np.mean(tractions) contact_points = np.where(1.0e-10 >= true_gap) contact_area = (np.size(contact_points)/2)/float(grid_size*grid_size) print('The contact area computed with the gap map is', contact_area) hertz_contact_size = (3 * pressure / (4 * curvature * effective_modulus))**(1. / 3.) print('Exact Hertz contact radius is', hertz_contact_size) hertz_area = np.pi * hertz_contact_size**2 print('Exact Hertz contact area is', hertz_area) area_error = np.abs(hertz_area - contact_area) / hertz_area print("Area error:", area_error) # Testing maximum pressure max_pressure = tractions.max() hertz_max_pressure = (6 * pressure * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi pressure_error = np.abs(hertz_max_pressure - max_pressure) \ / hertz_max_pressure print('Exact Hertz max pressure is', hertz_max_pressure) print('max pressure is', np.max(tractions)) print("Max pressure error:", pressure_error) # Testing displacements hertz_displacements = computeHertzDisplacement(effective_modulus, hertz_contact_size, hertz_max_pressure, grid_size) # Selecing only the points that are in contact contact_indexes = [(i, j) for i in range(grid_size) for j in range(grid_size) if true_gap[i, j] < 1e-12] # Displacements of bem are centered around the mean of the whole surface # and Hertz displacements are not centered, so we need to compute mean # on the contact zone for both arrays bem_mean = 0. hertz_mean = 0. for index in contact_indexes: bem_mean += true_displacements[index] hertz_mean += hertz_displacements[index] bem_mean /= len(contact_indexes) hertz_mean /= len(contact_indexes) # Correction applied when computing error correction = hertz_mean - bem_mean # Computation of error of displacement in contact zone error = 0. hertz_norm = 0. for index in contact_indexes: error += (hertz_displacements[index] - true_displacements[index] - correction)**2 hertz_norm += (hertz_displacements[index] - hertz_mean)**2 displacement_error = np.sqrt(error / hertz_norm) print("Displacement error (in contact zone): {}".format(displacement_error)) assert area_error < 1e-2 assert pressure_error < 1e-2 assert displacement_error < 2e-3 if __name__ == "__main__": test_hertz_disp() diff --git a/tests/test_tangential.py b/tests/test_tangential.py index 2a9f26e..ed986e0 100644 --- a/tests/test_tangential.py +++ b/tests/test_tangential.py @@ -1,102 +1,102 @@ # -*- coding: utf-8 -*- # @file # @section LICENSE # # Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # 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 . from __future__ import division, print_function import numpy as np from numpy.linalg import norm import tamaas as tm def test_pressure(tamaas_fixture): E = 1.0 nu = 0.5 mu = 0.5 n = 81 L = 1.0 p_target = np.array([0.5e-4, 0, 2e-4], dtype=tm.dtype) # g_target = np.array([-0.000223, 0, 9.99911], dtype=tm.dtype) x = np.linspace(-L/2, L/2, n, dtype=tm.dtype) y = np.copy(x) x, y = np.meshgrid(x, y, indexing="ij") r_sqrd = (x**2 + y**2) # Spherical contact R = 10.0 surface = R * np.sqrt((r_sqrd < R**2) * (1 - r_sqrd / R**2)) # Creating model - model = tm.ModelFactory.createModel(tm.model_type_surface_2d, + model = tm.ModelFactory.createModel(tm.model_type.surface_2d, [L, L], [n, n]) model.E, model.nu = E, nu p = model['traction'] # Theoretical solution Estar = E / (1.0 - nu**2) Fn = p_target[2] * L**2 Fx = p_target[0] * L**2 d_theory = np.power(3 / 4 * Fn / Estar / np.sqrt(R), 2/3) a_theory = np.sqrt(R * d_theory) c_theory = a_theory * (1 - Fx / (mu * Fn)) ** (1/3) p0_theory = np.power(6 * Fn * Estar**2 / np.pi**3 / R**2, 1/3) t1_theory = mu * p0_theory t2_theory = t1_theory * c_theory / a_theory # p_theory = p0_theory * np.sqrt((r_sqrd < a_theory**2) # * (1 - r_sqrd / a_theory**2)) t_theory = t1_theory * np.sqrt((r_sqrd < a_theory**2) * (1 - r_sqrd / a_theory**2)) t_theory = t_theory - t2_theory * np.sqrt((r_sqrd < c_theory**2) * (1 - r_sqrd / c_theory**2)) def assert_error(err): error = norm(p[..., 0] - t_theory) / norm(t_theory) print(error) assert error < err # Test Kato solver solver = tm.Kato(model, surface, 1e-12, mu) solver.max_iter = 1000 solver.solve(p_target, 100) assert_error(4e-2) # solver.solveRelaxed(g_target) # assert_error(1e-2) # # Test BeckTeboulle solver # solver = tm.BeckTeboulle(model, surface, 1e-12, mu) # solver.solve(g_target) # assert_error(1e-2) # Test Condat solver solver = tm.Condat(model, surface, 1e-12, mu) solver.max_iter = 5000 # or 10000 solver.solve(p_target) assert_error(7e-2) # 4e-2 for 10000 iterations # Test tangential Polonsky Keer solver solver = tm.PolonskyKeerTan(model, surface, 1e-12, mu) solver.max_iter = 1000 solver.solve(p_target) assert_error(4e-2) if __name__ == "__main__": test_pressure()