Page MenuHomec4science
Files F92819842

No OneTemporary
Actions

File Metadata

Created
Sat, Nov 23, 23:11

View Options

diff --git a/src/SConscript b/src/SConscript
index 59cb0b8..b4e6fec 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,114 +1,115 @@
import os
Import('main_env')
def prepend(path, list):
return [os.path.join(path, x) for x in list]
env = main_env
# Core
core_list = """
fft_plan_manager.cpp
fftransform.cpp
fftransform_fftw.cpp
grid.cpp
grid_hermitian.cpp
statistics.cpp
surface.cpp
tamaas.cpp
types.cpp
loop.cpp
""".split()
core_list = prepend('core', core_list)
core_list.append('tamaas_info.cpp')
# Lib roughcontact
generator_list = """
surface_generator.cpp
surface_generator_filter.cpp
surface_generator_filter_fft.cpp
isopowerlaw.cpp
""".split()
generator_list = prepend('surface', generator_list)
# Lib PERCOLATION
percolation_list = """
flood_fill.cpp
""".split()
percolation_list = prepend('percolation', percolation_list)
# BEM PERCOLATION
bem_list = """
bem_kato.cpp
bem_polonski.cpp
bem_gigi.cpp
bem_gigipol.cpp
bem_penalty.cpp
bem_uzawa.cpp
bem_fft_base.cpp
bem_functional.cpp
bem_meta_functional.cpp
elastic_energy_functional.cpp
exponential_adhesion_functional.cpp
squared_exponential_adhesion_functional.cpp
maugis_adhesion_functional.cpp
complimentary_term_functional.cpp
bem_grid.cpp
bem_grid_polonski.cpp
bem_grid_kato.cpp
bem_grid_teboulle.cpp
bem_grid_condat.cpp
""".split()
bem_list = prepend('bem', bem_list)
# Model
model_list = """
model.cpp
model_factory.cpp
model_type.cpp
model_template.cpp
westergaard.cpp
elastic_functional.cpp
meta_functional.cpp
adhesion_functional.cpp
elasto_plastic_model.cpp
ve_engine.cpp
kelvin.cpp
mindlin.cpp
""".split()
model_list = prepend('model', model_list)
# Solvers
solvers_list = """
contact_solver.cpp
polonsky_keer_rey.cpp
+kato.cpp
""".split()
solvers_list = prepend('solvers', solvers_list)
# GPU API
gpu_list = """
fftransform_cufft.cpp
""".split()
gpu_list = prepend('gpu', gpu_list)
# Assembling total list
rough_contact_list = \
core_list + generator_list + percolation_list + model_list + solvers_list + bem_list
# Adding GPU if needed
if env['backend'] == 'cuda':
rough_contact_list += gpu_list
# Generating dependency files
# env.AppendUnique(CXXFLAGS=['-MMD'])
# for file_name in rough_contact_list:
# obj = file_name.replace('.cpp', '.os')
# dep_file_name = file_name.replace('.cpp', '.d')
# env.SideEffect(dep_file_name, obj)
# env.ParseDepends(dep_file_name)
libTamaas = env.SharedLibrary('Tamaas', rough_contact_list)
Export('libTamaas')
diff --git a/src/solvers/kato.cpp b/src/solvers/kato.cpp
new file mode 100644
index 0000000..224650e
--- /dev/null
+++ b/src/solvers/kato.cpp
@@ -0,0 +1,312 @@
+/**
+ * @file
+ * @author Son Pham-Ba <son.phamba@epfl.ch>
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "kato.hh"
+#include "elastic_functional.hh"
+#include "loop.hh"
+#include "model_type.hh"
+#include "fft_plan_manager.hh"
+#include <iomanip>
+#include <iterator>
+/* -------------------------------------------------------------------------- */
+
+__BEGIN_TAMAAS__
+
+Kato::Kato(Model& model, const GridBase<Real>& surface,
+ Real tolerance, type variable_type,
+ type constraint_type)
+ : ContactSolver(model, surface, tolerance), variable_type(variable_type),
+ constraint_type(constraint_type), functional(model.getBEEngine()) {
+ if (model.getType() != model_type::basic_1d &&
+ model.getType() != model_type::basic_2d) {
+ TAMAAS_EXCEPTION(
+ "Model type is not compatible with Kato solver");
+ }
+
+ /// Gap is locally allocated
+ auto gap_ =
+ allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
+ model.getTraction().getNbComponents());
+ auto pressure_ = &model.getTraction();
+
+ switch (variable_type) {
+ case gap: {
+ primal = gap_;
+ dual = pressure_;
+ functional.addFunctionalTerm(new functional::ElasticFunctionalGap(
+ model.getBEEngine(), this->surface),
+ true);
+ break;
+ }
+ case pressure: {
+ primal = pressure_;
+ dual = gap_;
+ functional.addFunctionalTerm(new functional::ElasticFunctionalPressure(
+ model.getBEEngine(), this->surface),
+ true);
+ break;
+ }
+ }
+
+ search_direction =
+ allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
+ model.getTraction().getNbComponents());
+
+ projected_search_direction =
+ allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
+ model.getTraction().getNbComponents());
+}
+
+/* -------------------------------------------------------------------------- */
+
+Real Kato::solve(Real target) {
+ Real G = 0, Gold = 1, error = 0, error_norm = 1;
+ UInt n = 0;
+ bool delta = false;
+ *search_direction = 0;
+ //*dual = 0;
+
+ // Printing column headers
+ std::cout << std::setw(5) << "Iter"
+ << " " << std::setw(15) << "Cost_f"
+ << " " << std::setw(15) << "Error" << '\n'
+ << std::fixed;
+
+ // Setting uniform value if constraint
+ if (constraint_type == variable_type && std::abs(primal->sum()) <= 0)
+ *primal = target;
+ else if (constraint_type == variable_type)
+ *primal *= target / primal->mean();
+ else if (constraint_type != variable_type)
+ *primal = this->surface_stddev;
+
+ do {
+ // Computing functional gradient
+ functional.computeGradF(*primal, *dual);
+
+ 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);
+ } while (error > this->tolerance && n++ < this->max_iterations);
+
+ // Final update of dual variable
+ functional.computeGradF(*primal, *dual);
+
+ Real dual_min = dual->min();
+ *dual -= dual_min;
+
+ // Primal is pressure: need to make sure gap is positive
+ if (variable_type == pressure) {
+ model.getDisplacement() = *dual;
+ model.getDisplacement() += this->surface;
+ }
+ // Primal is gap: need to make sure dual is positive (pressure + adhesion)
+ else {
+ model.getDisplacement() = *primal;
+ model.getDisplacement() += this->surface;
+ this->model.solveDirichlet();
+ model.getTraction() -= dual_min;
+ }
+
+ return error;
+}
+
+/* -------------------------------------------------------------------------- */
+
+/**
+ * Computes \f$ \frac{1}{\mathrm{card}(\{p > 0\})} \sum_{\{p > 0\}}{f_i} \f$
+ */
+Real Kato::meanOnUnsaturated(const GridBase<Real>& field) const {
+ // Sum on unsaturated
+ Real sum = Loop::reduce<operation::plus>(
+ [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f : 0; },
+ *primal, field);
+ // Number of unsaturated points
+ UInt n_unsat = Loop::reduce<operation::plus>(
+ [] CUDA_LAMBDA(Real & p) -> UInt { return (p > 0); }, *primal);
+
+ return sum / n_unsat;
+}
+
+/* -------------------------------------------------------------------------- */
+
+/**
+ * Computes \f$ \sum_{\{p > 0\}}{f_i^2} \f$
+ */
+Real Kato::computeSquaredNorm(const GridBase<Real>& field) const {
+ return Loop::reduce<operation::plus>(
+ [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f * f : 0; },
+ *primal, field);
+}
+
+/* -------------------------------------------------------------------------- */
+
+/**
+ * Computes \f$ \tau = \frac{ \sum_{\{p > 0\}}{q_i't_i} }{ \sum_{\{p > 0\}}{r_i'
+ * t_i} } \f$
+ */
+Real Kato::computeCriticalStep(Real target) {
+ switch (variable_type) {
+ case gap:
+ model.getBEEngine().solveDirichlet(*search_direction,
+ *projected_search_direction);
+ break;
+ case pressure:
+ model.getBEEngine().solveNeumann(*search_direction,
+ *projected_search_direction);
+ break;
+ }
+
+ Real rbar = meanOnUnsaturated(*projected_search_direction);
+ if (constraint_type == variable_type)
+ *projected_search_direction -= rbar;
+ else {
+ *projected_search_direction += 2 * target + rbar;
+ }
+
+ Real num = Loop::reduce<operation::plus>(
+ [] CUDA_LAMBDA(Real & p, Real & q, Real & t) {
+ return (p > 0) ? q * t : 0;
+ },
+ *primal, *dual, *search_direction);
+ Real denum = Loop::reduce<operation::plus>(
+ [] CUDA_LAMBDA(Real & p, Real & r, Real & t) {
+ return (p > 0) ? r * t : 0;
+ },
+ *primal, *projected_search_direction, *search_direction);
+
+ return num / denum;
+}
+
+/* -------------------------------------------------------------------------- */
+
+/**
+ * Update steps:
+ * 1. \f$\mathbf{p} = \mathbf{p} - \tau \mathbf{t} \f$
+ * 2. Truncate all \f$p\f$ negative
+ * 3. For all points in \f$I_\mathrm{na} = \{p = 0 \land q < 0 \}\f$ do \f$p_i =
+ * p_i - \tau q_i\f$
+ */
+bool Kato::updatePrimal(Real step) {
+ UInt na_num = Loop::reduce<operation::plus>(
+ [step] CUDA_LAMBDA(Real & p, Real & q, Real & t) -> UInt {
+ p -= step * t; // Updating primal
+ if (p < 0)
+ p = 0; // Truncating neg values
+ if (p == 0 && q < 0) { // If non-admissible state
+ p -= step * q;
+ return 1;
+ } else
+ return 0;
+ },
+ *primal, *dual, *search_direction);
+
+ return na_num == 0;
+}
+
+/* -------------------------------------------------------------------------- */
+
+/**
+ * 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;
+}
+
+/* -------------------------------------------------------------------------- */
+
+/**
+ * Do \f$\mathbf{t} = \mathbf{q}' + \delta \frac{R}{R_\mathrm{old}}\mathbf{t}
+ * \f$
+ */
+void Kato::updateSearchDirection(Real factor) {
+ Loop::loop(
+ [factor] CUDA_LAMBDA(Real & p, Real & q, Real & t) {
+ t = (p > 0) ? q + factor * t : 0;
+ },
+ *primal, *dual, *search_direction);
+}
+
+/* -------------------------------------------------------------------------- */
+Kato::~Kato() {
+ switch (variable_type) {
+ case gap:
+ delete primal;
+ break;
+ case pressure:
+ delete dual;
+ break;
+ }
+
+ delete search_direction;
+ delete projected_search_direction;
+ FFTPlanManager::get().clean(); // TODO find a smarter way to deal with plans
+}
+
+/* -------------------------------------------------------------------------- */
+void Kato::addFunctionalTerm(functional::Functional * func) {
+ functional.addFunctionalTerm(func, false);
+}
+
+__END_TAMAAS__
+/* -------------------------------------------------------------------------- */
diff --git a/src/solvers/kato.hh b/src/solvers/kato.hh
new file mode 100644
index 0000000..bc7196d
--- /dev/null
+++ b/src/solvers/kato.hh
@@ -0,0 +1,77 @@
+/**
+ * @file
+ *
+ * @author Son Pham-Ba <son.phamba@epfl.ch>
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#ifndef __KATO__
+#define __KATO__
+/* -------------------------------------------------------------------------- */
+#include "contact_solver.hh"
+#include "meta_functional.hh"
+#include "tamaas.hh"
+/* -------------------------------------------------------------------------- */
+__BEGIN_TAMAAS__
+
+class Kato : public ContactSolver {
+public:
+ /// Types of algorithm (primal/dual) or constraint
+ enum type { gap, pressure };
+
+public:
+ /// Constructor
+ Kato(Model& model, const GridBase<Real>& surface, Real tolerance,
+ type variable_type, type constraint_type);
+ /// Destructor
+ ~Kato() override;
+
+public:
+ /// Solve
+ Real solve(Real target);
+ /// Mean on unsaturated constraint zone
+ Real meanOnUnsaturated(const GridBase<Real>& field) const;
+ /// Compute squared norm
+ Real computeSquaredNorm(const GridBase<Real>& field) const;
+ /// Update search direction
+ void updateSearchDirection(Real factor);
+ /// Compute critical step
+ Real computeCriticalStep(Real target = 0);
+ /// Update primal and check non-admissible state
+ bool updatePrimal(Real step);
+ /// Compute error/stopping criterion
+ Real computeError();
+ /// Add term to functional (fuck swig)
+ void addFunctionalTerm(::tamaas::functional::Functional * func);
+
+protected:
+ type variable_type, constraint_type;
+ GridBase<Real>* primal = nullptr;
+ GridBase<Real>* dual = nullptr;
+ GridBase<Real>* search_direction = nullptr;
+ GridBase<Real>* projected_search_direction = nullptr;
+ functional::MetaFunctional functional;
+};
+
+__END_TAMAAS__
+
+#endif // __KATO__

Event Timeline

c4science · Help