Page MenuHomec4science
Files F60100288

solvers.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Apr 27, 12:13

solvers.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2023 Lucas Frérot
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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 <pybind11/iostream.h>
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
/* -------------------------------------------------------------------------- */
using namespace py::literals;
class PyEPSolver : public EPSolver {
public:
using EPSolver::EPSolver;
// NOLINTNEXTLINE(readability-else-after-return)
void solve() override { PYBIND11_OVERLOAD_PURE(void, EPSolver, solve); }
void updateState() override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD(void, EPSolver, updateState);
}
};
class PyContactSolver : public ContactSolver {
public:
using ContactSolver::ContactSolver;
Real solve(std::vector<Real> load) override {
PYBIND11_OVERLOAD(Real, ContactSolver, solve, load);
}
Real solve(Real load) override {
PYBIND11_OVERLOAD(Real, ContactSolver, solve, load);
}
};
/* -------------------------------------------------------------------------- */
void wrapSolvers(py::module& mod) {
py::class_<ContactSolver, PyContactSolver>(mod, "ContactSolver")
.def(py::init<Model&, const GridBase<Real>&, Real>(),
py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
.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("tolerance", &ContactSolver::getTolerance,
&ContactSolver::setTolerance, "Solver tolerance")
.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_property_readonly("functional", &ContactSolver::getFunctional)
.def("addFunctionalTerm", &ContactSolver::addFunctionalTerm,
"Add a term to the contact functional to minimize")
.def("solve", py::overload_cast<std::vector<Real>>(&ContactSolver::solve),
"target_force"_a,
py::call_guard<py::scoped_ostream_redirect,
py::scoped_estream_redirect>(),
"Solve the contact for a mean traction/gap vector")
.def("solve", py::overload_cast<Real>(&ContactSolver::solve),
"target_normal_pressure"_a,
py::call_guard<py::scoped_ostream_redirect,
py::scoped_estream_redirect>(),
"Solve the contact for a mean normal pressure/gap");
py::class_<PolonskyKeerRey, ContactSolver> pkr(
mod, "PolonskyKeerRey",
"Main solver class for normal elastic contact problems. Its functional "
"can be customized to add an adhesion term, and its primal variable can "
"be set to either the gap or the pressure.");
// Need to export enum values before defining PKR constructor
py::enum_<PolonskyKeerRey::type>(pkr, "type")
.value("gap", PolonskyKeerRey::gap)
.value("pressure", PolonskyKeerRey::pressure)
.export_values();
pkr.def(py::init<Model&, const GridBase<Real>&, 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_<KatoSaturated, PolonskyKeerRey>(
mod, "KatoSaturated",
"Solver for pseudo-plasticity problems where the normal pressure is "
"constrained above by a saturation pressure \"pmax\"")
.def(py::init<Model&, const GridBase<Real>&, 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,
"Saturation normal pressure");
py::class_<Kato, ContactSolver> kato(mod, "Kato");
kato.def(py::init<Model&, const GridBase<Real>&, 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_<BeckTeboulle, ContactSolver> bt(mod, "BeckTeboulle");
bt.def(py::init<Model&, const GridBase<Real>&, 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_<Condat, ContactSolver> cd(
mod, "Condat",
"Main solver for frictional contact problems. It has no restraint on the "
"material properties or friction coefficient values, but solves an "
"associated version of the Coulomb friction law, which differs from the "
"traditional Coulomb friction in that the normal and tangential slip "
"components are coupled.");
cd.def(py::init<Model&, const GridBase<Real>&, 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_<PolonskyKeerTan, ContactSolver> pkt(mod, "PolonskyKeerTan");
pkt.def(py::init<Model&, const GridBase<Real>&, 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_<ToleranceManager>(
mod, "_tolerance_manager",
"Manager object for the tolereance of nonlinear plasticity solvers. "
"Decreases the solver tolerance by geometric progression.")
.def(py::init<Real, Real, Real>(), "start_tol"_a, "end_tol"_a, "rate"_a);
py::class_<EPSolver, PyEPSolver>(
mod, "EPSolver", "Mother class for nonlinear plasticity solvers")
.def(py::init<Residual&>(), "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_<DFSANESolver, EPSolver>(mod, "_DFSANESolver")
.def(py::init<Residual&>(), "residual"_a, py::keep_alive<1, 2>())
.def(py::init([](Residual& res, Model&) {
TAMAAS_DEPRECATE("Solver(residual, model)", "Solver(residual)");
return std::make_unique<DFSANESolver>(res);
}),
"residual"_a, "model"_a, py::keep_alive<1, 2>());
py::class_<EPICSolver>(
mod, "EPICSolver",
"Main solver class for elastic-plastic contact problems")
.def(py::init<ContactSolver&, EPSolver&, Real, Real>(),
"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<Real>{pressure});
},
"normal_pressure"_a,
py::call_guard<py::scoped_ostream_redirect,
py::scoped_estream_redirect>(),
"Solves the EP contact with a relaxed fixed-point scheme. Adjust "
"the relaxation parameter to help convergence.")
.def(
"acceleratedSolve",
[](EPICSolver& solver, Real pressure) {
return solver.acceleratedSolve(std::vector<Real>{pressure});
},
"normal_pressure"_a,
py::call_guard<py::scoped_ostream_redirect,
py::scoped_estream_redirect>(),
"Solves the EP contact with an accelerated fixed-point scheme. May "
"not converge!");
}
} // namespace wrap
} // namespace tamaas

Event Timeline

c4science · Help