Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/solvers/kato.cpp
356b5974ef35features/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,
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 = this->_gap.get();
pressure = &model.getTraction();
this->functional.addFunctionalTerm(new functional::ElasticFunctionalPressure(
model.getBEEngine(), this->surface),
true);
}
/* -------------------------------------------------------------------------- */
Real Kato::solve(GridBase<Real>& p0) {
Real error = 0, error_norm = 1, cost_f = 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);
auto test = pressure->getInternalData();
test[2] = 1.0;
// do {
*gap = 0;
functional.computeGradF(*pressure, *gap);
*pressure -= *gap;
enforcePressureConstraints(p0);
// 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;
}
/* -------------------------------------------------------------------------- */
/**
* Projects $\vec{p}$ on $\mathcal{C}$ and $\mathcal{D}$.
*/
void Kato::enforcePressureConstraints(GridBase<Real>& p0) {
UInt n = 10;
for (UInt i = 0; i < n; i++) {
enforcePressureMean(p0);
enforcePressureCoulomb();
}
}
/* -------------------------------------------------------------------------- */
void Kato::enforcePressureMean(GridBase<Real>& p0) {
Grid<Real, 1> mean({pressure->getNbComponents()}, 1);
// initialize mean with same size as p0
computeMean(*pressure, mean);
auto begin_it(pressure->begin());
auto end_it(pressure->end());
#pragma omp parallel for
for (auto it = begin_it; it < end_it; ++it) {
UInt i = it - begin_it;
*it += p0(i % 3)
- mean(i % 3);
}
}
/* -------------------------------------------------------------------------- */
void Kato::enforcePressureCoulomb() {
#pragma omp parallel for
for (auto i = pressure->begin(); i < pressure->end(); i += 3) {
auto it = i;
auto p1 = *it;
auto p2 = *(++it);
auto p_n = *(++it);
auto p_t_sqrd = p1 * p1 + p2 * p2;
}
}
/* -------------------------------------------------------------------------- */
/**
* Compute mean of the field taking each component separately.
*/
void Kato::computeMean(GridBase<Real>& field, GridBase<Real>& 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__
/* -------------------------------------------------------------------------- */
c4science · Help