Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/solvers/kato.cpp
44877ef96c21maxwell_viscoelastic…

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();//&model.getDisplacement();
pressure = &model.getTraction();
this->functional.addFunctionalTerm(new functional::ElasticFunctionalPressure(
model.getBEEngine(), this->surface),
true);
}
/* -------------------------------------------------------------------------- */
Real Kato::solve(GridBase<Real>& p0) {
Real cost = 0;
UInt n = 0;
Real beta = 0;
Grid<Real, 1> lambda({pressure->getNbPoints()}, 1);
// 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);
*gap += surface;
// *gap -= surface; // already computed in computeGradF
// *pressure -= *gap;
// enforcePressureConstraints(p0);
// beta = computeBeta();
// computeLambda(lambda, beta);
// cost = computeCost(lambda);
// printState(n, cost, cost);
// } while (cost > this->tolerance && n++ < this->max_iterations);
// computeFinalGap(beta);
model.getDisplacement() = *gap;
return cost;
}
/* -------------------------------------------------------------------------- */
/**
* Projects $\vec{p}$ on $\mathcal{C}$ and $\mathcal{D}$.
*/
void Kato::enforcePressureConstraints(GridBase<Real>& p0) {
UInt n = 100;
for (UInt i = 0; i < n; i++) {
enforcePressureMean(p0);
enforcePressureCoulomb();
}
}
/* -------------------------------------------------------------------------- */
void Kato::enforcePressureMean(GridBase<Real>& p0) {
Grid<Real, 1> corr({pressure->getNbComponents()}, 1);
// initialize mean with same size as p0
computeMean(*pressure, corr);
corr -= p0;
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 -= corr(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;
if (p_n >= 0 && p_t_sqrd <= mu * mu * p_n * p_n) {
continue;
} else if (p_n <= 0 && p_t_sqrd <= p_n * p_n / mu / mu) {
it = i;
*it = 0;
*(++it) = 0;
*(++it) = 0;
} else {
auto p_t = std::sqrt(p_t_sqrd);
Real k = (p_n + mu * p_t) / (1 + mu * mu);
it = i;
*it = k * mu * p1 / p_t;
*(++it) = k * mu * p2 / p_t;
*(++it) = k;
}
}
}
/* -------------------------------------------------------------------------- */
/**
* 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();
}
/* -------------------------------------------------------------------------- */
Real Kato::computeBeta() {
Real beta = -std::numeric_limits<Real>::infinity();
#pragma omp parallel for
for (auto i = gap->begin(); i < gap->end(); i += 3) {
auto it = i;
auto g1 = *it;
auto g2 = *(++it);
auto g_n = *(++it);
auto g_t = std::sqrt(g1 * g1 + g2 * g2);
auto beta_temp = mu * g_t - g_n;
beta = (beta_temp > beta) ? beta_temp : beta;
}
return beta;
}
/* -------------------------------------------------------------------------- */
void Kato::computeLambda(GridBase<Real>& lambda, Real beta) {
#pragma omp parallel for
for (auto i = gap->begin(); i < gap->end(); i += 3) {
auto it = i;
auto g1 = *it;
auto g2 = *(++it);
auto g_n = *(++it);
auto g_t = std::sqrt(g1 * g1 + g2 * g2);
lambda((i - gap->begin()) / 3) = g_n - mu * g_t + beta;
}
}
/* -------------------------------------------------------------------------- */
Real Kato::computeCost(GridBase<Real>& lambda) {
UInt N = pressure->getNbPoints();
Grid<Real, 1> P_n({N}, 1);
Grid<Real, 1> P_c({N}, 1);
Grid<Real, 1> eta({N}, 1);
#pragma omp parallel for
for (UInt i = 0; i < N; ++i) {
auto p1 = (*pressure)(3 * i);
auto p2 = (*pressure)(3 * i + 1);
auto p_n = (*pressure)(3 * i + 2);
auto p_t = std::sqrt(p1 * p1 + p2 * p2);
auto g1 = (*gap)(3 * i);
auto g2 = (*gap)(3 * i + 1);
auto g_t = std::sqrt(g1 * g1 + g2 * g2);
P_n(i) = p_n;
P_c(i) = mu * p_n - p_t;
eta(i) = g_t;
}
return P_n.dot(lambda) + P_c.dot(eta);
}
/* -------------------------------------------------------------------------- */
void Kato::computeFinalGap(Real beta) {
#pragma omp parallel for
for (UInt i = 0; i < gap->getNbPoints(); ++i) {
(*gap)(3 * i + 2) += surface(3 * i + 2) + beta;
}
}
/* -------------------------------------------------------------------------- */
void Kato::addFunctionalTerm(functional::Functional * func) {
functional.addFunctionalTerm(func, false);
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
c4science · Help