diff --git a/python/wrap/core.cpp b/python/wrap/core.cpp index a587d36..2b83ff9 100644 --- a/python/wrap/core.cpp +++ b/python/wrap/core.cpp @@ -1,109 +1,112 @@ /** * @file * LICENSE * * Copyright (©) 2016-2021 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 "logger.hh" #include "statistics.hh" #include "wrap.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ namespace wrap { using namespace py::literals; /* -------------------------------------------------------------------------- */ template void wrapStatistics(py::module& mod) { auto name = makeDimensionName("Statistics", dim); py::class_>(mod, name.c_str()) .def_static("computePowerSpectrum", &Statistics::computePowerSpectrum, - py::return_value_policy::move) - .def_static("computeAutocorrelation", - &Statistics::computeAutocorrelation, - py::return_value_policy::move) - .def_static("computeMoments", &Statistics::computeMoments) + py::return_value_policy::move, "Compute PSD of surface") + .def_static( + "computeAutocorrelation", &Statistics::computeAutocorrelation, + py::return_value_policy::move, "Compute autocorrelation of surface") + .def_static("computeMoments", &Statistics::computeMoments, + "Compute spectral moments") .def_static("computeSpectralRMSSlope", - &Statistics::computeSpectralRMSSlope) - .def_static("computeRMSHeights", &Statistics::computeRMSHeights) + &Statistics::computeSpectralRMSSlope, + "Compute hrms' in Fourier space") + .def_static("computeRMSHeights", &Statistics::computeRMSHeights, + "Compute hrms") .def_static("contact", &Statistics::contact, "tractions"_a, "perimeter"_a = 0, "Compute the (corrected) contact area. Permieter is the " "total contact perimeter in number of segments."); } /* -------------------------------------------------------------------------- */ void wrapCore(py::module& mod) { // Exposing logging facility py::enum_(mod, "LogLevel") .value("debug", LogLevel::debug) .value("info", LogLevel::info) .value("warning", LogLevel::warning) .value("error", LogLevel::error); mod.def("set_log_level", [](LogLevel level) { Logger::setLevel(level); }); py::class_(mod, "Logger") .def(py::init<>()) .def("get", [](Logger& logger, LogLevel level) -> Logger& { logger.get(level); return logger; }) .def("__lshift__", [](Logger& logger, std::string msg) -> Logger& { auto level = logger.getWishLevel(); if (level >= Logger::getCurrentLevel() and not(mpi::rank() != 0 and level == LogLevel::info)) { py::print(msg, "file"_a = py::module::import("sys").attr("stderr")); } return logger; }); wrapStatistics<1>(mod); wrapStatistics<2>(mod); mod.def("to_voigt", [](const Grid& field) { - Logger().get(LogLevel::warning) - << "tamaas.to_voigt deprecated. Use tamaas.compute.to_voigt\n"; + TAMAAS_DEPRECATE("tamaas.to_voigt()", "tamaas.compute.to_voigt()"); + if (field.getNbComponents() == 9) { Grid voigt(field.sizes(), 6); Loop::loop( [] CUDA_LAMBDA(MatrixProxy in, SymMatrixProxy out) { out.symmetrize(in); }, range>(field), range>(voigt)); return voigt; } else TAMAAS_EXCEPTION("Wrong number of components to symmetrize"); }, "Convert a 3D tensor field to voigt notation", py::return_value_policy::copy); } } // namespace wrap /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp index 87fb854..4f160c7 100644 --- a/python/wrap/model.cpp +++ b/python/wrap/model.cpp @@ -1,431 +1,451 @@ /** * @file * LICENSE * * Copyright (©) 2016-2021 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) - .def("setParameters", - [](functional::AdhesionFunctional& f, - const std::map& m) { - TAMAAS_DEPRECATE("setParameters()", "the parameters property"); - f.setParameters(m); - }); + .def("setParameters", [](functional::AdhesionFunctional& f, + const std::map& m) { + 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") .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); }) .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); + 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_1·N_2·N_3) complexity, may cause float " + "overflow/underflow") + .value("cutoff", integration_method::cutoff, + "Approximation, O(sqrt(N_1^2+N_2^2)·N_3^2) 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", - [](py::object) { return trait::dimension; }) + [](py::object) { return trait::dimension; }, + "Dimension of computational domain") .def_property_readonly_static( - "components", [](py::object) { return trait::components; }) + "components", [](py::object) { return trait::components; }, + "Number of components of vector fields") .def_property_readonly_static( "boundary_dimension", - [](py::object) { return trait::boundary_dimension; }) - .def_property_readonly_static("voigt", - [](py::object) { return trait::voigt; }) + [](py::object) { return trait::boundary_dimension; }, + "Dimension of boundary of computational domain") + .def_property_readonly_static( + "voigt", [](py::object) { return trait::voigt; }, + "Number of components of symmetrical tensor fields") .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); + .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") .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) + .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) - .def("solveDirichlet", &Model::solveDirichlet) - .def("dump", &Model::dump) - .def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>()) + .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(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", [](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("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[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 " "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/solvers.cpp b/python/wrap/solvers.cpp index 42ca273..eddb3ff 100644 --- a/python/wrap/solvers.cpp +++ b/python/wrap/solvers.cpp @@ -1,188 +1,189 @@ /** * @file * LICENSE * * Copyright (©) 2016-2021 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") .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("addFunctionalTerm", &ContactSolver::addFunctionalTerm, + "Add a term to the contact functional to minimize") .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/src/core/surface_complex.hh b/src/core/surface_complex.hh deleted file mode 100644 index c9f136b..0000000 --- a/src/core/surface_complex.hh +++ /dev/null @@ -1,159 +0,0 @@ -/** - * @file - * LICENSE - * - * Copyright (©) 2016-2021 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 SURFACE_COMPLEX_H -#define SURFACE_COMPLEX_H -/* -------------------------------------------------------------------------- */ -#include "grid_hermitian.hh" -#include "surface.hh" -#include - -namespace tamaas { - -/** - * @deprecated should use GridHermitian instead - * @brief Square 2D complex surface - */ -template -class SurfaceComplex : public GridHermitian { - -public: - SurfaceComplex(UInt a, Real L); - using GridHermitian::GridHermitian; - using GridHermitian::operator=; - using GridHermitian::operator(); - - //! make a surface real by taking norm as real part and imaginary 0 - void makeItRealBySquare(); - //! make it real by summing the imaginary and real parts - void makeItRealBySum(); - //! make a surface real by taking norm as real part and imaginary 0 - void makeItRealByAbs(); - //! get real part - Surface real() const; - //! get imaginary part - Surface imag() const; - - UInt size() const { return this->n[0]; } - UInt getL() const { return 1.; } - - void setGridSize(UInt s); - - virtual const SurfaceComplex& getWrappedNumpy() { return *this; }; - -#define SURF_OPERATOR(op) \ - inline void operator op(const Surface& other); \ - inline void operator op(const SurfaceComplex& other) - - SURF_OPERATOR(+=); - SURF_OPERATOR(*=); - SURF_OPERATOR(-=); - SURF_OPERATOR(/=); -#undef SURF_OPERATOR -}; - -/* -------------------------------------------------------------------------- */ - -template -void SurfaceComplex::setGridSize(UInt s) { - this->n[0] = s; - this->n[1] = s / 2 + 1; - this->resize(this->n); -} - -/* -------------------------------------------------------------------------- */ - -template -Surface SurfaceComplex::real() const { - Surface res(this->size(), this->getL()); - UInt n = this->size(); - for (UInt i = 0; i < n; ++i) - for (UInt j = 0; j < n; ++j) { - res(i, j) = (*this)(i, j).real(); - } - return res; -} -/* -------------------------------------------------------------------------- */ - -template -Surface SurfaceComplex::imag() const { - Surface res(this->size(), this->getL()); - UInt n = this->size(); - for (UInt i = 0; i < n; ++i) - for (UInt j = 0; j < n; ++j) { - res(i, j) = (*this)(i, j).imag(); - } - return res; -} - -/* -------------------------------------------------------------------------- */ -template -void SurfaceComplex::makeItRealBySquare() { - Loop::loop([] CUDA_LAMBDA(complex & x) { x *= thrust::conj(x); }, *this); -} - -/* -------------------------------------------------------------------------- */ -template -void SurfaceComplex::makeItRealByAbs() { - Loop::loop([] CUDA_LAMBDA(complex & x) { x = thrust::abs(x); }, *this); -} - -/* -------------------------------------------------------------------------- */ -template -void SurfaceComplex::makeItRealBySum() { - Loop::loop([] CUDA_LAMBDA(complex & x) { x = x.real() + x.imag(); }, - *this); -} - -/* -------------------------------------------------------------------------- */ -template -SurfaceComplex::SurfaceComplex(UInt a, Real /*L*/) : GridHermitian() { - this->n = GridHermitian::hermitianDimensions(std::array{a, a}); - this->resize(this->n); -} - -/* -------------------------------------------------------------------------- */ -#define SURF_OP_IMPL(op) \ - template \ - inline void SurfaceComplex::operator op(const Surface& other) { \ - TAMAAS_ASSERT(this->n[0] == other.sizes()[0] && \ - this->n[1] == other.sizes()[1] / 2 + 1, \ - "surface size does not match"); \ - _Pragma("omp parallel for") for (UInt i = 0; i < this->n[0]; \ - i++) for (UInt j = 0; j < this->n[1]; \ - j++) (*this)(i, j)op other(i, \ - j); \ - } \ - template \ - inline void SurfaceComplex::operator op(const SurfaceComplex& other) { \ - Grid, 2>::operator op(other); \ - } - -SURF_OP_IMPL(+=) -SURF_OP_IMPL(*=) -SURF_OP_IMPL(-=) -SURF_OP_IMPL(/=) -#undef SURF_OP_IMPL - -} // namespace tamaas - -#endif /* SURFACE_COMPLEX_H */ diff --git a/src/model/kelvin.cpp b/src/model/kelvin.cpp index ee230f4..a615a12 100644 --- a/src/model/kelvin.cpp +++ b/src/model/kelvin.cpp @@ -1,143 +1,140 @@ /** * @file * LICENSE * * Copyright (©) 2016-2021 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 "kelvin.hh" #include "logger.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* Constructors */ /* -------------------------------------------------------------------------- */ template Kelvin::Kelvin(Model* model) : VolumePotential(model) { setIntegrationMethod(integration_method::linear, 1e-12); } template void Kelvin::setIntegrationMethod(integration_method method, Real cutoff) { this->method = method; this->cutoff = cutoff; Logger logger; if (this->method == integration_method::linear) { logger.get(LogLevel::debug) << TAMAAS_DEBUG_MSG("Setting linear integration method"); this->initialize(dtrait::template source_components, dtrait::template out_components, this->model->getDiscretization()[0]); } else { logger.get(LogLevel::debug) << TAMAAS_DEBUG_MSG( "Setting cutoff integration method (cutoff " << this->cutoff << ')'); this->initialize(dtrait::template source_components, dtrait::template out_components, 1); } auto max_q = Loop::reduce( [] CUDA_LAMBDA(VectorProxy qv) { return qv.l2norm(); }, range>( this->wavevectors)); if (this->method == integration_method::linear and not std::isfinite(std::exp(max_q * this->model->getSystemSize()[0]))) logger.get(LogLevel::warning) << "Probable overflow of integral computation (consider " "changing integration method to integration_method::cutoff or " "compiling with real_type='long double')\n"; } /* -------------------------------------------------------------------------- */ /* Operator implementation */ /* -------------------------------------------------------------------------- */ template void Kelvin::applyIf(GridBase& source, GridBase& out, filter_t pred) const { KelvinInfluence kelvin(this->model->getShearModulus(), this->model->getPoissonRatio()); parent::transformSource(source, pred); // Reset buffer values for (auto&& layer : this->out_buffer) layer = 0; if (method == integration_method::linear) { linearIntegral(out, kelvin); } else { cutoffIntegral(out, kelvin); } } template void Kelvin::linearIntegral(GridBase& out, KelvinInfluence& kelvin) const { detail::KelvinHelper helper; helper.applyIntegral(this->source_buffer, this->out_buffer, this->wavevectors, this->model->getSystemSize().front(), kelvin); helper.applyFreeTerm(this->source_buffer, this->out_buffer, kelvin); helper.makeFundamentalGreatAgain(this->out_buffer); parent::transformOutput( [](auto&& out_buffer, auto layer) -> typename parent::BufferType& { return out_buffer[layer]; }, out); } template void Kelvin::cutoffIntegral(GridBase& out, KelvinInfluence& kelvin) const { detail::KelvinHelper helper; auto func = [&](auto&& out_buffer, auto layer) -> typename parent::BufferType& { auto&& out = out_buffer.front(); out = 0; helper.applyIntegral(this->source_buffer, out, layer, this->wavevectors, this->model->getSystemSize().front(), cutoff, kelvin); helper.applyFreeTerm(this->source_buffer[layer], out, kelvin); helper.makeFundamentalGreatAgain(out); return out; }; parent::transformOutput(func, out); } -/* -------------------------------------------------------------------------- - */ +/* -------------------------------------------------------------------------- */ /* Template instanciation */ -/* -------------------------------------------------------------------------- - */ +/* -------------------------------------------------------------------------- */ template class Kelvin; template class Kelvin; template class Kelvin; -/* -------------------------------------------------------------------------- - */ +/* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/src/surface/isopowerlaw.cpp b/src/surface/isopowerlaw.cpp index d8f89e5..08f54f6 100644 --- a/src/surface/isopowerlaw.cpp +++ b/src/surface/isopowerlaw.cpp @@ -1,101 +1,101 @@ /** * @file * LICENSE * * Copyright (©) 2016-2021 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 /* -------------------------------------------------------------------------- */ namespace tamaas { template void Isopowerlaw::computeFilter( GridHermitian& filter_coefficients) const { Filter::computeFilter( [this] CUDA_LAMBDA(Complex & coeff, VectorProxy q) { coeff = (*this)(q); }, filter_coefficients); } /* -------------------------------------------------------------------------- */ template Real Isopowerlaw::rmsHeights() const { return std::sqrt(moments()[0]); } /* -------------------------------------------------------------------------- */ /** * Analytical moments, cf. Yastrebov et al. (2015) * "From infinitesimal to full contact between rough surfaces: Evolution * of the contact area", appendix A */ template <> std::vector Isopowerlaw<2>::moments() const { std::map T; T[0] = 2 * M_PI; T[2] = M_PI; T[4] = 3 * M_PI / 4.; Real q1 = static_cast(this->q1); // shadowing this->q1 Real xi = q0 / q1; Real zeta = q2 / q1; Real A = std::pow(q1, 2 * (hurst + 1)); std::vector moments; moments.reserve(3); for (UInt q : {0, 2, 4}) { Real m = A * T[q] * std::pow(q1, q - 2. * hurst) * ((1 - std::pow(xi, q + 2)) / (q + 2.) + (std::pow(zeta, q - 2. * hurst) - 1) / (q - 2. * hurst)); moments.push_back(m); } return moments; } template <> std::vector Isopowerlaw<1>::moments() const { TAMAAS_EXCEPTION("Moments have not been implemented for 1D surfaces"); } /* -------------------------------------------------------------------------- */ template Real Isopowerlaw::alpha() const { - auto m = moments(); + const auto m = moments(); return m[0] * m[2] / (m[1] * m[1]); } /* -------------------------------------------------------------------------- */ template Real Isopowerlaw::rmsSlopes() const { // 2 pi because moments are computed with k instead of q // and slopes should be computed with q return 2 * M_PI * std::sqrt(2 * moments()[1]); } template class Isopowerlaw<1>; template class Isopowerlaw<2>; } // namespace tamaas