diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp index eb13ac1..e5771ea 100644 --- a/python/wrap/solvers.cpp +++ b/python/wrap/solvers.cpp @@ -1,73 +1,73 @@ /** * @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" /* -------------------------------------------------------------------------- */ __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>(), - "model"_a, "surface"_a, "tolerance"_a, + 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) .def("computeError", &Kato::computeError) .def("addFunctionalTerm", &Kato::addFunctionalTerm); } } // namespace wrap __END_TAMAAS__ diff --git a/src/solvers/kato.cpp b/src/solvers/kato.cpp index f5d4f32..674b406 100644 --- a/src/solvers/kato.cpp +++ b/src/solvers/kato.cpp @@ -1,139 +1,177 @@ /** * @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 "model_type.hh" #include "fft_plan_manager.hh" #include #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ -Kato::Kato(Model& model, const GridBase& surface, Real tolerance) - : ContactSolver(model, surface, tolerance), - functional(model.getBEEngine()) { +Kato::Kato(Model& model, const GridBase& surface, Real tolerance, + Real mu) + : ContactSolver(model, surface, tolerance), mu(mu) { if (model.getType() != model_type::surface_2d) { TAMAAS_EXCEPTION( "Model type is not compatible with Kato solver"); } - /// Gap is locally allocated - // gap = - // allocateGrid(model.getType(), model.getDiscretization(), - // model.getTraction().getNbComponents()); - gap = &model.getDisplacement(); + gap = this->_gap.get(); pressure = &model.getTraction(); - functional.addFunctionalTerm(new functional::ElasticFunctionalGap( - model.getBEEngine(), this->surface), - true); + this->functional.addFunctionalTerm(new functional::ElasticFunctionalPressure( + model.getBEEngine(), this->surface), + true); } /* -------------------------------------------------------------------------- */ Real Kato::solve(GridBase& p0) { - Real error = 0, error_norm = 1; + Real error = 0, error_norm = 1, cost_f = 0; UInt n = 0; - + Grid mean({pressure->getNbComponents()}, 1); // initialize mean with same size as p0 + // Printing column headers std::cout << std::setw(5) << "Iter" << " " << std::setw(15) << "Cost_f" << " " << std::setw(15) << "Error" << '\n' << std::fixed; - + // pressure->uniformSetComponents(p0); + auto test = pressure->getInternalData(); + test[512*256*3+256*3+2] = 1.0; // do { - // // Computing functional gradient - // functional.computeGradF(*primal, *dual); + *gap = 0; + functional.computeGradF(*pressure, *gap); + *pressure -= *gap; + computeMean(*pressure, mean); + enforcePressureConstraints(); // Real dbar = meanOnUnsaturated(*dual); // // Enforcing dual constraint via gradient // if (constraint_type != variable_type) { // *dual += 2 * target + dbar; // } else { // // Centering dual on primal > 0 // *dual -= dbar; // } // // Computing gradient norm // G = computeSquaredNorm(*dual); // // Updating search direction (conjugate gradient) // updateSearchDirection(delta * G / Gold); // Gold = G; // // Computing critical step // Real tau = computeCriticalStep(target); // // Update primal variable // delta = updatePrimal(tau); // // We should scale to constraint // if (constraint_type == variable_type) // *primal *= target / primal->mean(); // error = computeError() / error_norm; // Real cost_f = functional.computeF(*primal, *dual); - // printState(n, cost_f, error); + printState(n, cost_f, error); // } while (error > this->tolerance && n++ < this->max_iterations); - + return error; } /* -------------------------------------------------------------------------- */ +/** + * Projects $\vec{p}$ on $\mathcal{C}$ and $\mathcal{D}$. + */ +void Kato::enforcePressureConstraints() { + UInt n = 10; + + for (UInt i = 0; i < n; i++) { + + } +} + +/* -------------------------------------------------------------------------- */ + +/** + * Compute mean of the field taking each component separately. + */ +void Kato::computeMean(GridBase& field, GridBase& mean) { + if (!(mean.dataSize() == field.getNbComponents())) + TAMAAS_EXCEPTION("Size of mean variable is not correct."); + + mean = 0; + auto mean_data = mean.getInternalData(); + auto begin_it(field.begin()); + auto end_it(field.end()); + +#pragma omp parallel for + for (auto it = begin_it; it < end_it; ++it) { + UInt i = it - begin_it; + mean_data[i % field.getNbComponents()] += *it; + } + + mean /= field.getNbPoints(); +} + +/* -------------------------------------------------------------------------- */ + /** * Error is based on \f$ \sum{p_i q_i} \f$ */ Real Kato::computeError() { // /// Making sure dual respects constraint // *dual -= dual->min(); // Real error = primal->dot(*dual); // Real norm = 1; // if (variable_type == pressure) // norm = std::abs(primal->sum() * this->surface_stddev); // else // norm = std::abs(dual->sum() * this->surface_stddev); // norm *= primal->getNbPoints(); // return std::abs(error) / norm; return 0.0; } /* -------------------------------------------------------------------------- */ void Kato::addFunctionalTerm(functional::Functional * func) { functional.addFunctionalTerm(func, false); } __END_TAMAAS__ /* -------------------------------------------------------------------------- */ diff --git a/src/solvers/kato.hh b/src/solvers/kato.hh index 5622854..b76fb68 100644 --- a/src/solvers/kato.hh +++ b/src/solvers/kato.hh @@ -1,59 +1,59 @@ /** * @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__ #define __KATO__ /* -------------------------------------------------------------------------- */ #include "contact_solver.hh" #include "meta_functional.hh" #include "tamaas.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class Kato : public ContactSolver { public: /// Constructor - Kato(Model& model, const GridBase& surface, Real tolerance); - /// Destructor - ~Kato() override; + Kato(Model& model, const GridBase& surface, Real tolerance, Real mu); public: /// Solve Real solve(GridBase& p0); /// Compute error/stopping criterion + void enforcePressureConstraints(); + void computeMean(GridBase& field, GridBase& mean); Real computeError(); - /// Add term to functional (fuck swig) - void addFunctionalTerm(::tamaas::functional::Functional * func); + /// Add term to functional + void addFunctionalTerm(functional::Functional * func); protected: GridBase* gap = nullptr; GridBase* pressure = nullptr; - functional::MetaFunctional functional; + Real mu = 0; }; __END_TAMAAS__ #endif // __KATO__