diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp index 50d90d0..3aac3e7 100644 --- a/python/wrap/solvers.cpp +++ b/python/wrap/solvers.cpp @@ -1,98 +1,99 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "wrap.hh" #include "contact_solver.hh" #include "polonsky_keer_rey.hh" #include "kato.hh" #include "beck_teboulle.hh" #include "condat.hh" #include "polonsky_keer_tan.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ namespace wrap { /* -------------------------------------------------------------------------- */ using namespace py::literals; /* -------------------------------------------------------------------------- */ void wrapSolvers(py::module& mod) { py::class_(mod, "ContactSolver") .def(py::init&, Real>(), "model"_a, "surface"_a, "tolerance"_a) .def("setMaxIterations", &ContactSolver::setMaxIterations, "max_iter"_a) .def("setDumpFrequency", &ContactSolver::setDumpFrequency, "dump_freq"_a) .def("addFunctionalTerm", &ContactSolver::addFunctionalTerm); py::class_ pkr(mod, "PolonskyKeerRey"); pkr.def(py::init&, Real, PolonskyKeerRey::type, PolonskyKeerRey::type>(), "model"_a, "surface"_a, "tolerance"_a, "primal_type"_a, "constraint_type"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &PolonskyKeerRey::solve, "target"_a) .def("computeError", &PolonskyKeerRey::computeError); py::enum_(pkr, "type") .value("gap", PolonskyKeerRey::gap) .value("pressure", PolonskyKeerRey::pressure) .export_values(); 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); + .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("computeCost", &PolonskyKeerTan::computeCost); + .def("solveTresca", &PolonskyKeerTan::solveTresca, "p0"_a) + .def("computeCost", &PolonskyKeerTan::computeCost, "use_tresca"_a=false); } } // namespace wrap __END_TAMAAS__ diff --git a/src/solvers/kato.cpp b/src/solvers/kato.cpp index b6cee8c..fbd7c34 100644 --- a/src/solvers/kato.cpp +++ b/src/solvers/kato.cpp @@ -1,539 +1,659 @@ /** * @file * * @author Son Pham-Ba * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "kato.hh" #include "elastic_functional.hh" #include "loop.hh" #include "fft_plan_manager.hh" #include #include +#include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ Kato::Kato(Model& model, const GridBase& surface, Real tolerance, Real mu) : ContactSolver(model, surface, tolerance), engine(model.getBEEngine()), mu(mu) { if (model.getType() != model_type::surface_1d && model.getType() != model_type::surface_2d) { TAMAAS_EXCEPTION( "Model type is not compatible with Kato solver"); } gap = this->_gap.get(); // locally allocated pressure = &model.getTraction(); N = pressure->getNbPoints(); if (model.getType() == model_type::surface_1d) { initSurfaceWithComponents(); } else { initSurfaceWithComponents(); } } /* -------------------------------------------------------------------------- */ Real Kato::solve(GridBase& p0, UInt proj_iter) { if (p0.getNbPoints() != pressure->getNbComponents()) { TAMAAS_EXCEPTION( "Target mean pressure does not have the right number of components"); } Real cost = 0; switch (model.getType()) { case model_type::surface_1d: cost = solveTmpl(p0, proj_iter); break; case model_type::surface_2d: cost = solveTmpl(p0, proj_iter); break; default: break; } return cost; } template Real Kato::solveTmpl(GridBase& p0, UInt proj_iter) { constexpr UInt comp = model_type_traits::components; Real cost = 0; UInt n = 0; // Printing column headers std::cout << std::setw(5) << "Iter" << " " << std::setw(15) << "Cost_f" << " " << std::setw(15) << "Error" << '\n' << std::fixed; pressure->uniformSetComponents(p0); do { computeGradient(); *pressure -= *gap; enforcePressureConstraints(p0, proj_iter); cost = computeCost(); printState(n, cost, cost); } while (cost > this->tolerance && n++ < this->max_iterations); computeFinalGap(); return cost; } /* -------------------------------------------------------------------------- */ Real Kato::solveRelaxed(GridBase& g0) { if (g0.getNbPoints() != pressure->getNbComponents()) { TAMAAS_EXCEPTION( "Target mean gap does not have the right number of components"); } Real cost = 0; switch (model.getType()) { case model_type::surface_1d: cost = solveRelaxedTmpl(g0); break; case model_type::surface_2d: cost = solveRelaxedTmpl(g0); break; default: break; } return cost; } template Real Kato::solveRelaxedTmpl(GridBase& g0) { constexpr UInt comp = model_type_traits::components; Real cost = 0; UInt n = 0; // Printing column headers std::cout << std::setw(5) << "Iter" << " " << std::setw(15) << "Cost_f" << " " << std::setw(15) << "Error" << '\n' << std::fixed; *pressure = 0; do { engine.solveNeumann(*pressure, *gap); addUniform(*gap, g0); *gap -= *surfaceComp; *pressure -= *gap; enforcePressureCoulomb(); cost = computeCost(); printState(n, cost, cost); } while (cost > this->tolerance && n++ < this->max_iterations); computeFinalGap(); return cost; } /* -------------------------------------------------------------------------- */ Real Kato::solveRegularized(GridBase& p0, Real r) { if (p0.getNbPoints() != pressure->getNbComponents()) { TAMAAS_EXCEPTION( "Target mean pressure does not have the right number of components"); } Real cost = 0; switch (model.getType()) { case model_type::surface_1d: cost = solveRegularizedTmpl(p0, r); break; case model_type::surface_2d: cost = solveRegularizedTmpl(p0, r); break; default: break; } return cost; } template Real Kato::solveRegularizedTmpl(GridBase& p0, Real r) { constexpr UInt comp = model_type_traits::components; Real cost = 0; UInt n = 0; // Printing column headers std::cout << std::setw(5) << "Iter" << " " << std::setw(15) << "Cost_f" << " " << std::setw(15) << "Error" << '\n' << std::fixed; pressure->uniformSetComponents(p0); do { // enforcePressureMean(p0); engine.solveNeumann(*pressure, *gap); *gap -= *surfaceComp; // Impose zero tangential displacement in non-sliding zone UInt count_static = 0; Vector g_static = Loop::stridedReduce( [&] CUDA_LAMBDA(VectorProxy&& g, VectorProxy&& p) -> Vector { VectorProxy p_T = p(0); Real p_T_norm = p_T.l2norm(); Real p_N = p(comp - 1); if (0.99 * mu * p_N > p_T_norm) { // non-sliding contact count_static ++; return g; // to compute mean of g_T } else { return 0; } }, *gap, *pressure); g_static /= count_static != 0 ? count_static : 1; g_static(comp - 1) = 0; Loop::stridedLoop( [this, r, g_static] CUDA_LAMBDA(VectorProxy&& p, VectorProxy&& g) { // Add frictional term to gradient of functional g -= g_static; - Vector _g = g; // copy + // Vector _g = g; // copy VectorProxy g_T (g(0)); VectorProxy g_N (g(comp - 1)); Real g_T_norm = g_T.l2norm(); // g_N += mu * regularize(g_T_norm, r) * g_T_norm; g_N += mu * g_T_norm; // Update pressure with gradient // _g *= 0.1; p -= g; // Truncate negative normal pressure VectorProxy p_T (p(0)); VectorProxy p_N (p(comp - 1)); if (p_N(0) < 0) p_N = 0; // Set tangential pressure p_T = g_T; if (g_T_norm != 0) p_T *= - mu * p_N(0) * regularize(g_T_norm, r) / g_T_norm; }, *pressure, *gap); // enforcePressureMean(p0); // enforcePressureCoulomb(); enforcePressureConstraints(p0, 50); cost = computeCost(); printState(n, cost, cost); } while (std::abs(cost) > this->tolerance && n++ < this->max_iterations); computeFinalGap(); return cost; } /* -------------------------------------------------------------------------- */ template void Kato::initSurfaceWithComponents() { constexpr UInt comp = model_type_traits::components; surfaceComp = allocateGrid(type, model.getDiscretization(), comp); *surfaceComp = 0; Loop::stridedLoop( [] CUDA_LAMBDA(Real& s, VectorProxy&& sc) { sc(comp - 1) = s; }, surface, *surfaceComp); } /* -------------------------------------------------------------------------- */ template -void Kato::computeGradient() { +void Kato::computeGradient(bool use_tresca) { engine.solveNeumann(*pressure, *gap); *gap -= *surfaceComp; // Impose zero tangential displacement in non-sliding zone UInt count_static = 0; - Vector g_static = Loop::stridedReduce( - [&] CUDA_LAMBDA(VectorProxy&& g, - VectorProxy&& p) -> Vector { - VectorProxy p_T = p(0); - Real p_T_norm = p_T.l2norm(); - Real p_N = p(comp - 1); + Vector g_static; - if (0.99 * mu * p_N > p_T_norm) { // non-sliding contact - count_static ++; - return g; // to compute mean of g_T - } else { - return 0; - } - }, - *gap, *pressure); + if (!use_tresca) { // with Coulomb friction + + g_static = Loop::stridedReduce( + [&] CUDA_LAMBDA(VectorProxy&& g, + VectorProxy&& p) -> Vector { + VectorProxy p_T = p(0); + Real p_T_norm = p_T.l2norm(); + Real p_N = p(comp - 1); + + if (0.99 * mu * p_N > p_T_norm) { // non-sliding contact + count_static ++; + return g; // to compute mean of g_T + } else { + return 0; + } + }, + *gap, *pressure); + + } else { // with Tresca friction + + g_static = Loop::stridedReduce( + [&] CUDA_LAMBDA(VectorProxy&& g, + VectorProxy&& p) -> Vector { + VectorProxy p_T = p(0); + Real p_N = p(comp - 1); + Real p_T_norm = p_T.l2norm(); + + if (0.99 * mu > p_T_norm && p_N > 0) { // non-sliding contact + count_static ++; + return g; // to compute mean of g_T + } else { + return 0; + } + }, + *gap, *pressure); + + } if (count_static != 0) { g_static /= count_static; // g_static(comp - 1) = 0; } else { // if no sticking zone, mean computed on sliding zone count_static = 0; Real g_N_static = Loop::stridedReduce( [&] CUDA_LAMBDA(VectorProxy&& g, VectorProxy&& p) { if (p(comp - 1) > 0) { count_static ++; return g(comp - 1); } else { return 0.0; } }, *gap, *pressure); g_static(comp - 1) = g_N_static / count_static; } - // Add frictionnal term to functional - Loop::stridedLoop( - [this, g_static] CUDA_LAMBDA(VectorProxy&& g) { - g -= g_static; - VectorProxy g_T = g(0); - Real g_T_norm = g_T.l2norm(); - g(comp - 1) += mu * g_T_norm; // Frictionnal work - }, - *gap); + if (!use_tresca) { // with Coulomb friction + + Loop::stridedLoop( + [this, g_static] CUDA_LAMBDA(VectorProxy&& g) { + g -= g_static; + VectorProxy g_T = g(0); + Real g_T_norm = g_T.l2norm(); + g(comp - 1) += mu * g_T_norm; // Frictionnal work + }, + *gap); + + } else { // with Tresca friction + *gap -= g_static; + + // Loop::stridedLoop( + // [this, g_static] CUDA_LAMBDA(VectorProxy&& g, + // VectorProxy&& p) { + // g -= g_static; + // VectorProxy g_T = g(0); + // Real g_T_norm = g_T.l2norm(); + // Real p_N = p(comp - 1); + // if (p_N > 0) + // g(comp - 1) += mu/p_N * g_T_norm; // Frictionnal work + // }, + // *gap, *pressure); + + } // Loop::stridedLoop( // [this] CUDA_LAMBDA(VectorProxy&& g, VectorProxy&& p) { // VectorProxy g_T = g(0); // Real g_T_norm = g_T.l2norm(); // VectorProxy g_N = g(comp - 1); // VectorProxy p_T = p(0); // Real p_T_norm = p_T.l2norm(); // VectorProxy p_N = p(comp - 1); // if (p_T_norm != 0) { // g_T = p_T; // g_T *= -g_T_norm / p_T_norm; // } // g_N += mu * g_T_norm; // Frictionnal work // }, // *gap, *pressure); } /* -------------------------------------------------------------------------- */ /** * Projects $\vec{p}$ on $\mathcal{C}$ and $\mathcal{D}$. */ template void Kato::enforcePressureConstraints(GridBase& p0, UInt proj_iter) { for (UInt i = 0; i < proj_iter; i++) { enforcePressureMean(p0); enforcePressureCoulomb(); } } /* -------------------------------------------------------------------------- */ template void Kato::enforcePressureMean(GridBase& p0) { Vector corr = computeMean(*pressure); VectorProxy _p0 = p0(0); corr -= _p0; *pressure -= corr; } /* -------------------------------------------------------------------------- */ template void Kato::enforcePressureCoulomb() { Loop::stridedLoop( [this] CUDA_LAMBDA(VectorProxy&& p) { VectorProxy p_T(p(0)); Real p_N = p(comp - 1); Real p_T_sqrd= p_T.l2squared(); // Projection normale au cône de friction bool cond1 = (p_N >= 0 && p_T_sqrd <= mu * mu * p_N * p_N); bool cond2 = (p_N <= 0 && p_T_sqrd <= p_N * p_N / mu / mu); if (cond2) { p_T = 0; p(comp - 1) = 0; } else if (!cond1) { Real p_T_norm = std::sqrt(p_T_sqrd); Real k = (p_N + mu * p_T_norm) / (1 + mu * mu); p_T *= k * mu / p_T_norm; p(comp - 1) = k; } }, *pressure); } /* -------------------------------------------------------------------------- */ +template +void Kato::enforcePressureTresca() { + Loop::stridedLoop( + [this] CUDA_LAMBDA(VectorProxy&& p) { + VectorProxy p_T(p(0)); + Real p_N = p(comp - 1); + Real p_T_norm = p_T.l2norm(); + + if (p_N < 0) { + p_T = 0; + p(comp - 1) = 0; + } else if (p_T_norm > mu) { // TODO: replace name of variable mu + p_T *= mu / p_T_norm; + // p(comp - 1) += p_T_norm - mu; + } + }, + *pressure); +} + +template void Kato::enforcePressureTresca<2>(); +template void Kato::enforcePressureTresca<3>(); + +/* -------------------------------------------------------------------------- */ + /** * Compute mean of the field taking each component separately. */ template Vector Kato::computeMean(GridBase& field) { Vector mean = Loop::stridedReduce( [] CUDA_LAMBDA(VectorProxy&& f) -> Vector { return f; }, field); mean /= N; return mean; } /* -------------------------------------------------------------------------- */ template void Kato::addUniform(GridBase& field, GridBase& vec) { VectorProxy _vec(vec(0)); field += _vec; } /* -------------------------------------------------------------------------- */ -Real Kato::computeCost() { +Real Kato::computeCost(bool use_tresca) { UInt N = pressure->getNbPoints(); Real beta = 0; Grid lambda({N}, 1); Grid eta({N}, 1); Grid p_N({N}, 1); Grid p_C({N}, 1); switch (model.getType()) { case model_type::surface_1d: - beta = computeBeta(); - computeValuesForCost(beta, lambda, eta, p_N, p_C); + if (!use_tresca) { + beta = computeBeta(); + computeValuesForCost(beta, lambda, eta, p_N, p_C); + } else { + computeValuesForCostTresca(lambda, eta, p_N, p_C); + } break; + case model_type::surface_2d: - beta = computeBeta(); - computeValuesForCost(beta, lambda, eta, p_N, p_C); + if (!use_tresca) { + beta = computeBeta(); + computeValuesForCost(beta, lambda, eta, p_N, p_C); + } else { + computeValuesForCostTresca(lambda, eta, p_N, p_C); + } break; + default: break; } return p_N.dot(lambda) + p_C.dot(eta); } /* -------------------------------------------------------------------------- */ template Real Kato::computeBeta() { constexpr UInt comp = model_type_traits::components; return Loop::stridedReduce( [this] CUDA_LAMBDA(VectorProxy&& g) { VectorProxy g_T(g(0)); Real g_N = g(comp - 1); Real g_T_norm = g_T.l2norm(); return mu * g_T_norm - g_N; }, *gap); } /* -------------------------------------------------------------------------- */ template void Kato::computeValuesForCost(Real beta, GridBase& lambda, GridBase& eta, GridBase& p_N, GridBase& p_C) { constexpr UInt comp = model_type_traits::components; Loop::stridedLoop( [this, beta] CUDA_LAMBDA(VectorProxy&& p, VectorProxy&& g, Real& lambda_, Real& eta_, Real& p_N_, Real& p_C_) { VectorProxy g_T(g(0)); Real g_N = g(comp - 1); Real g_T_norm = g_T.l2norm(); - lambda_ = g_N - mu * g_T_norm + beta; + lambda_ = g_N - mu * g_T_norm + beta; eta_ = g_T_norm; VectorProxy p_T(p(0)); Real p_N = p(comp - 1); Real p_T_norm = p_T.l2norm(); p_N_ = p(comp - 1); p_C_ = mu * p_N - p_T_norm; }, *pressure, *gap, lambda, eta, p_N, p_C); } /* -------------------------------------------------------------------------- */ +template +void Kato::computeValuesForCostTresca(GridBase& lambda, GridBase& eta, + GridBase& p_N, GridBase& p_C) { + constexpr UInt comp = model_type_traits::components; + + Real g_N_min = Loop::stridedReduce( + [this] CUDA_LAMBDA(VectorProxy&& g) { + return g(comp - 1); + }, + *gap); + + Real p_C_min = Loop::stridedReduce( + [this] CUDA_LAMBDA(VectorProxy&& p) { + VectorProxy p_T(p(0)); + Real p_T_norm = p_T.l2norm(); + return mu - p_T_norm; + }, + *pressure); + + Loop::stridedLoop( + [this, g_N_min, p_C_min] CUDA_LAMBDA(VectorProxy&& p, + VectorProxy&& g, + Real& lambda_, + Real& eta_, + Real& p_N_, + Real& p_C_) { + VectorProxy g_T(g(0)); + Real g_N = g(comp - 1); + Real g_T_norm = g_T.l2norm(); + lambda_ = g_N - g_N_min; + eta_ = g_T_norm; + + VectorProxy p_T(p(0)); + Real p_N = p(comp - 1); + Real p_T_norm = p_T.l2norm(); + p_N_ = p_N; + p_C_ = (p_N > 0) ? mu - p_T_norm - p_C_min : 0; + }, + *pressure, *gap, lambda, eta, p_N, p_C); +} + +/* -------------------------------------------------------------------------- */ + template void Kato::computeFinalGap() { engine.solveNeumann(*pressure, *gap); *gap -= *surfaceComp; Real g_N_min = Loop::stridedReduce( [] CUDA_LAMBDA(VectorProxy&& g) { return g(comp - 1); }, *gap); Grid g_shift({comp}, 1); g_shift = 0; g_shift(comp - 1) = -g_N_min; *gap += *surfaceComp; addUniform(*gap, g_shift); model.getDisplacement() = *gap; } /* -------------------------------------------------------------------------- */ Real Kato::regularize(Real x, Real r) { Real xr = x / r; return xr / (1 + std::abs(xr)); } __END_TAMAAS__ /* -------------------------------------------------------------------------- */ diff --git a/src/solvers/kato.hh b/src/solvers/kato.hh index e488b08..9a88b32 100644 --- a/src/solvers/kato.hh +++ b/src/solvers/kato.hh @@ -1,112 +1,119 @@ /** * @file * * @author Son Pham-Ba * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __KATO_HH__ #define __KATO_HH__ /* -------------------------------------------------------------------------- */ #include "contact_solver.hh" #include "meta_functional.hh" #include "model_type.hh" #include "static_types.hh" #include "tamaas.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class Kato : public ContactSolver { public: /// Constructor Kato(Model& model, const GridBase& surface, Real tolerance, Real mu); public: /// Solve Real solve(GridBase& p0, UInt proj_iter); /// Solve relaxed problem Real solveRelaxed(GridBase& g0); /// Solve regularized problem Real solveRegularized(GridBase& p0, Real r); /// Compute cost function - Real computeCost(); + Real computeCost(bool use_tresca = false); private: /// Template for solve function template Real solveTmpl(GridBase& p0, UInt proj_iter); /// Template for solveRelaxed function template Real solveRelaxedTmpl(GridBase& g0); /// Template for solveRegularized function template Real solveRegularizedTmpl(GridBase& p0, Real r); protected: /// Creates surfaceComp form surface template void initSurfaceWithComponents(); /// Compute gradient of functional template - void computeGradient(); + void computeGradient(bool use_tresca = false); /// Project pressure on friction cone template void enforcePressureConstraints(GridBase& p0, UInt proj_iter); /// Project on C template void enforcePressureMean(GridBase& p0); /// Project on D template void enforcePressureCoulomb(); + /// Project on D (Tresca) + template + void enforcePressureTresca(); /// Comupte mean value of field template Vector computeMean(GridBase& field); /// Add vector to each point of field template void addUniform(GridBase& field, GridBase& vec); /// Regularization function with factor r (0 -> unregugularized) Real regularize(Real x, Real r); - /// Compute mean gap + /// Compute shift template Real computeBeta(); /// Compute grids of dual and primal variables template void computeValuesForCost(Real beta, GridBase& lambda, GridBase& eta, GridBase& p_N, GridBase& p_C); + /// Compute dual and primal variables with Tresca friction + template + void computeValuesForCostTresca(GridBase& lambda, + GridBase& eta, GridBase& p_N, GridBase& p_C); /// Compute total displacement template void computeFinalGap(); protected: BEEngine& engine; GridBase* gap = nullptr; GridBase* pressure = nullptr; std::unique_ptr> surfaceComp = nullptr; Real mu = 0; UInt N = 0; // number of points }; __END_TAMAAS__ #endif // __KATO_HH__ diff --git a/src/solvers/polonsky_keer_tan.cpp b/src/solvers/polonsky_keer_tan.cpp index 5be7aa2..47c09ba 100644 --- a/src/solvers/polonsky_keer_tan.cpp +++ b/src/solvers/polonsky_keer_tan.cpp @@ -1,304 +1,332 @@ /** * @file * * @author Son Pham-Ba * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "polonsky_keer_tan.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ PolonskyKeerTan::PolonskyKeerTan(Model& model, const GridBase& surface, Real tolerance, Real mu) : Kato(model, surface, tolerance, mu) { search_direction = allocateGrid(model.getType(), model.getDiscretization(), model.getTraction().getNbComponents()); search_direction_backup = allocateGrid(model.getType(), model.getDiscretization(), model.getTraction().getNbComponents()); projected_search_direction = allocateGrid(model.getType(), model.getDiscretization(), model.getTraction().getNbComponents()); } /* -------------------------------------------------------------------------- */ Real PolonskyKeerTan::solve(GridBase& p0) { if (p0.getNbPoints() != pressure->getNbComponents()) { TAMAAS_EXCEPTION( "Target mean pressure does not have the right number of components"); } Real cost = 0; switch (model.getType()) { case model_type::surface_1d: cost = solveTmpl(p0); break; case model_type::surface_2d: cost = solveTmpl(p0); break; default: break; } return cost; } +/* -------------------------------------------------------------------------- */ + +Real PolonskyKeerTan::solveTresca(GridBase& p0) { + if (p0.getNbPoints() != pressure->getNbComponents()) { + TAMAAS_EXCEPTION( + "Target mean pressure does not have the right number of components"); + } + + Real cost = 0; + + switch (model.getType()) { + case model_type::surface_1d: + cost = solveTmpl(p0, true); + break; + case model_type::surface_2d: + cost = solveTmpl(p0, true); + break; + default: + break; + } + + return cost; +} + template -Real PolonskyKeerTan::solveTmpl(GridBase& p0) { +Real PolonskyKeerTan::solveTmpl(GridBase& p0, bool use_tresca) { // Printing column headers std::cout << std::setw(5) << "Iter" << " " << std::setw(15) << "Cost_f" << " " << std::setw(15) << "Error" << '\n' << std::fixed; constexpr UInt comp = model_type_traits::components; Real cost = 0.0; UInt n = 0; pressure->uniformSetComponents(p0); Real R_old = 1.0; *search_direction = 0.0; do { // Enforce external condition (should be at the end) enforcePressureMean(p0); // Compute functional gradient - computeGradient(); + computeGradient(use_tresca); // Compute search direction Real R = computeSquaredNorm(*gap); *search_direction *= R / R_old; *search_direction += *gap; R_old = R; // Compute step size Real tau = computeStepSize(false); // Update pressure *search_direction *= tau; *pressure -= *search_direction; // Enforce constraints - enforcePressureCoulomb(); + if (!use_tresca) { + enforcePressureCoulomb(); + } else { + enforcePressureTresca(); + } - cost = computeCost(); + cost = computeCost(use_tresca); printState(n, cost, cost); } while (cost > this->tolerance && n++ < this->max_iterations); computeFinalGap(); return cost; } // Original algorithm // template // Real PolonskyKeerTan::solveTmpl(GridBase& p0) { // // Printing column headers // std::cout << std::setw(5) << "Iter" // << " " << std::setw(15) << "Cost_f" // << " " << std::setw(15) << "Error" << '\n' // << std::fixed; // constexpr UInt comp = model_type_traits::components; // Real cost = 0.0; // UInt n = 0; // pressure->uniformSetComponents(p0); // Real R_old = 1.0; // Real delta = 0.0; // *search_direction = 0.0; // do { // // Enforce external condition (should be at the end) // enforcePressureMean(p0); // // Compute functional gradient // computeGradient(); // Vector gap_mean = computeMean(*gap, true); // for (UInt i = 0; i < comp - 1; i++) gap_mean(i) = 0; // *gap -= gap_mean; // // Compute search direction // Real R = computeSquaredNorm(*gap); // *search_direction *= delta * R / R_old; // *search_direction += *gap; // R_old = R; // truncateSearchDirection(true); // // Compute step size // Real tau = computeStepSize(true); // // Update pressure // *search_direction *= tau; // *pressure -= *search_direction; // // Enforce constraints // enforcePressureCoulomb(); // // Empty set of inadmissible gaps // UInt na_count = Loop::stridedReduce( // [tau] CUDA_LAMBDA(VectorProxy&& p, // VectorProxy&& g) { // if (p(comp - 1) == 0.0 && g(comp - 1) < 0.0) { // Vector _g = g; // _g *= tau; // p -= _g; // return 1; // } else { // return 0; // } // }, // *pressure, *gap); // delta = (na_count > 0) ? 0.0 : 1.0; // // Enforce external condition // // enforcePressureMean(p0); // cost = computeCost(); // printState(n, cost, cost); // } while (cost > this->tolerance && n++ < this->max_iterations); // computeFinalGap(); // return cost; // } /* -------------------------------------------------------------------------- */ template void PolonskyKeerTan::enforcePressureMean(GridBase& p0) { Vector pressure_mean = computeMean(*pressure, false); // *pressure -= pressure_mean; // addUniform(*pressure, p0); // for (UInt i = 0; i < comp; i++) // if (pressure_mean(i) == 0) pressure_mean(i) = 1.0; // *pressure /= pressure_mean; // VectorProxy _p0(p0(0)); // *pressure *= _p0; Loop::stridedLoop( [pressure_mean, p0] CUDA_LAMBDA(VectorProxy&& p) { for (UInt i = 0; i < comp - 1; i++) p(i) += p0(i) - pressure_mean(i); p(comp - 1) *= p0(comp - 1) / pressure_mean(comp - 1); }, *pressure); // Loop::stridedLoop( // [pressure_mean, p0] CUDA_LAMBDA(VectorProxy&& p) { // for (UInt i = 0; i < comp - 1; i++) // p(i) *= p0(i) / (pressure_mean(i) != 0 ? pressure_mean(i) : 1); // p(comp - 1) += p0(comp - 1) - pressure_mean(comp - 1); // }, // *pressure); } /* -------------------------------------------------------------------------- */ template Vector PolonskyKeerTan::computeMean(GridBase& field, bool on_c) { UInt count = 0; Vector mean = Loop::stridedReduce( [&] CUDA_LAMBDA(VectorProxy&& f, VectorProxy&& p) -> Vector { if ((!on_c) || p(comp - 1) > 0.0) { count ++; return f; } else { return 0; } }, field, *pressure); mean /= count; return mean; } /* -------------------------------------------------------------------------- */ Real PolonskyKeerTan::computeSquaredNorm(GridBase& field) { return Loop::reduce( [] CUDA_LAMBDA(Real& f) { return f * f; }, field); } /* -------------------------------------------------------------------------- */ template void PolonskyKeerTan::truncateSearchDirection(bool on_c) { if (!on_c) return; Loop::stridedLoop( [] CUDA_LAMBDA(VectorProxy&& t, VectorProxy&& p) { if (p(comp - 1) == 0.0) t = 0.0; }, *search_direction, *pressure); } /* -------------------------------------------------------------------------- */ template Real PolonskyKeerTan::computeStepSize(bool on_c) { engine.solveNeumann(*search_direction, *projected_search_direction); Vector r_mean = computeMean(*projected_search_direction, on_c); *projected_search_direction -= r_mean; Real num = Loop::stridedReduce( [] CUDA_LAMBDA(VectorProxy&& q, VectorProxy&& t) { Real dot = 0.0; for (UInt i = 0; i < comp; i++) dot += q(i) * t(i); return dot; }, *gap, *search_direction); Real denum = Loop::stridedReduce( [] CUDA_LAMBDA(VectorProxy&& r, VectorProxy&& t) { Real dot = 0.0; for (UInt i = 0; i < comp; i++) dot += r(i) * t(i); return dot; }, *projected_search_direction, *search_direction); return num / denum; } __END_TAMAAS__ /* -------------------------------------------------------------------------- */ \ No newline at end of file diff --git a/src/solvers/polonsky_keer_tan.hh b/src/solvers/polonsky_keer_tan.hh index bb003c8..1f31ab4 100644 --- a/src/solvers/polonsky_keer_tan.hh +++ b/src/solvers/polonsky_keer_tan.hh @@ -1,71 +1,73 @@ /** * @file * * @author Son Pham-Ba * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __POLONSKY_KEER_TAN_HH__ #define __POLONSKY_KEER_TAN_HH__ /* -------------------------------------------------------------------------- */ #include "kato.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class PolonskyKeerTan : public Kato { public: /// Constructor PolonskyKeerTan(Model& model, const GridBase& surface, Real tolerance, Real mu); public: - /// Solve + /// Solve with Coulomb friction Real solve(GridBase& p0); + /// Solve with Tresca friction + Real solveTresca(GridBase& p0); private: /// Template for solve function template - Real solveTmpl(GridBase& p0); + Real solveTmpl(GridBase& p0, bool use_tresca = false); /// Enforce pressure mean template void enforcePressureMean(GridBase& p0); /// Compute mean of field (only on I_c) template Vector computeMean(GridBase& field, bool on_c); /// Compute squared norm Real computeSquaredNorm(GridBase& field); /// Restrict search direction on I_c template void truncateSearchDirection(bool on_c); /// Compute optimal step size (only on I_c) template Real computeStepSize(bool on_c); private: std::unique_ptr> search_direction = nullptr; std::unique_ptr> search_direction_backup = nullptr; std::unique_ptr> projected_search_direction = nullptr; }; __END_TAMAAS__ #endif // __POLONSKY_KEER_TAN_HH__