Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/solvers/kato.cpp
561cbfbcee77features/declarative

Last Change
Raw File

kato.cpp
View Options

/**
* @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)
: ContactSolver(model, surface, tolerance),
functional(model.getBEEngine()) {
if (model.getType() != model_type::surface_2d) {
TAMAAS_EXCEPTION(
"Model type is not compatible with Kato solver");
}
/// Gap is locally allocated
// gap =
// allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
// model.getTraction().getNbComponents());
gap = &model.getDisplacement();
pressure = &model.getTraction();
functional.addFunctionalTerm(new functional::ElasticFunctionalGap(
model.getBEEngine(), this->surface),
true);
}
/* -------------------------------------------------------------------------- */
Real Kato::solve(GridBase<Real>& p0) {
Real error = 0, error_norm = 1;
UInt n = 0;
// Printing column headers
std::cout << std::setw(5) << "Iter"
<< " " << std::setw(15) << "Cost_f"
<< " " << std::setw(15) << "Error" << '\n'
<< std::fixed;
// 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);
return error;
}
/* -------------------------------------------------------------------------- */
/**
* 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__
/* -------------------------------------------------------------------------- */
c4science · Help