Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/bem/bem_grid_kato.cpp
83b12befa682maxwell_viscoelastic…

Last Change
Raw File

bem_grid_kato.cpp
View Options

/**
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016 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 "bem_grid_kato.hh"
#include "surface_statistics.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
BemGridKato::BemGridKato(Surface<Real>& surface)
: BemGrid(surface), lambda(surface.size(), surface.getL()), beta(0) {
this->max_iter = 10000;
}
/* -------------------------------------------------------------------------- */
BemGridKato::~BemGridKato() {}
/* -------------------------------------------------------------------------- */
void BemGridKato::computeTrueTractions() {
// nothing to do here, tractions are our main unknown
}
void BemGridKato::computeTrueDisplacements() {
// nothing to do here, tractions are our main unknown
if (!this->true_displacements)
this->true_displacements = new Grid<Real, 2>(
this->displacements.sizes(), this->displacements.getNbComponents());
*this->true_displacements = this->displacements;
const UInt N = this->true_displacements->getNbPoints();
Real *d_start = this->true_displacements->getInternalData(),
*g_start = this->gaps.getInternalData();
legacy::VectorProxy<Real> d(nullptr, 3), g(nullptr, 3);
#pragma omp parallel for firstprivate(d, g)
for (UInt i = 0; i < N; i++) {
d.setPointer(d_start + i * d.dataSize());
g.setPointer(g_start + i * g.dataSize());
d(0) = g(0);
d(1) = g(1);
d(2) = g(2) + this->surface(i) + beta;
}
}
/* -------------------------------------------------------------------------- */
Real BemGridKato::computeFunctional() {
Real *p_start = this->tractions.getInternalData(),
*u_start = this->displacements.getInternalData();
legacy::VectorProxy<Real> p(nullptr, 3), u(nullptr, 3);
Real E = 0;
const UInt N = this->tractions.getNbPoints();
#pragma omp parallel for reduction(+ : E) firstprivate(p, u)
for (UInt i = 0; i < N; i++) {
p.setPointer(p_start + i * p.dataSize());
u.setPointer(u_start + i * u.dataSize());
E += 0.5 * p.dot(u) - p(2) * this->surface(i);
}
return E;
}
/* -------------------------------------------------------------------------- */
Real BemGridKato::computeLambda() {
Real* g_start = this->gaps.getInternalData();
legacy::VectorProxy<Real> g(nullptr, 3), g_T(nullptr, 2);
const UInt N = this->gaps.getNbPoints();
const UInt nb_components = this->gaps.getNbComponents();
Real beta = std::numeric_limits<Real>::max();
// Computing lambda and beta at the same time
#pragma omp parallel for reduction(min : beta) firstprivate(g, g_T)
for (UInt i = 0; i < N; i++) {
g.setPointer(g_start + i * nb_components);
g_T.setPointer(g_start + i * nb_components);
lambda(i) = g(2) - mu * g_T.l2norm();
if (lambda(i) < beta)
beta = lambda(i);
}
return beta;
}
/* -------------------------------------------------------------------------- */
void BemGridKato::updateGradient(Grid<Real, 1>& q_vec) {
legacy::VectorProxy<Real> q(q_vec.getInternalData(), 3), g(nullptr, 3);
const UInt N = gaps.getNbPoints();
#pragma omp parallel for firstprivate(g)
for (UInt i = 0; i < N; i++) {
g.setPointer(gaps.getInternalData() + i * g.dataSize());
g += q;
}
}
/* -------------------------------------------------------------------------- */
void BemGridKato::computeEquilibrium(Real tol, const Grid<Real, 1>& g_target) {
this->computeInfluence();
Real f = 0;
UInt n_iter = 0;
TAMAAS_ASSERT(g_target.dataSize() == 3,
"Imposed gap must have only 3 components");
do {
this->computeDisplacementsFromTractions();
this->computeGaps(this->displacements);
updateGradient(const_cast<Grid<Real, 1>&>(g_target));
beta = computeLambda();
this->tractions -= this->gaps;
BemGrid::projectOnFrictionCone(this->tractions, mu);
lambda -= beta;
f = computeStoppingCriterion();
printFriendlyMessage(n_iter, f);
} while (f > tol && n_iter++ < this->max_iter);
}
/* -------------------------------------------------------------------------- */
/// Stopping criterion is othogonality between lambda and pressure
Real BemGridKato::computeStoppingCriterion() {
Real max = std::numeric_limits<Real>::min();
Real ortho = 0;
Real *p_start = this->tractions.getInternalData(),
*g_start = this->gaps.getInternalData();
legacy::VectorProxy<Real> p(nullptr, 3), p_T(nullptr, 2), g_T(nullptr, 2);
const UInt N = this->tractions.dataSize() / this->tractions.getNbComponents();
const UInt nb_components = this->tractions.getNbComponents();
const Real hrms = SurfaceStatistics::computeSpectralStdev(this->surface);
#pragma omp parallel firstprivate(p, g_T, p_T)
{
#pragma omp for reduction(max : max)
for (UInt i = 0; i < this->tractions.dataSize(); i++) {
if (max < this->tractions(i))
max = this->tractions(i);
}
#pragma omp for reduction(+ : ortho)
for (UInt i = 0; i < N; i++) {
p.setPointer(p_start + i * nb_components);
p_T.setPointer(p_start + i * nb_components);
g_T.setPointer(g_start + i * nb_components);
ortho += std::abs(p(2) * lambda(i)) +
std::abs((mu * p(2) - p_T.l2norm()) * g_T.l2norm());
}
} // end #pragma omp parallel
return ortho / (max * hrms);
}
/* -------------------------------------------------------------------------- */
void BemGridKato::printFriendlyMessage(UInt niter, Real conv_crit) {
if (niter % this->dump_freq == 0) {
UInt prec = std::cout.precision();
std::cout << niter << " " << std::scientific << conv_crit << " "
<< computeFunctional() << std::endl;
std::cout << std::fixed;
std::cout.precision(prec);
}
// if (std::isnan(conv_crit))
// std::cerr << "Something went wrong in the solve phase !" << std::endl;
}
} // namespace tamaas
c4science · Help