Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/bem/bem_grid_kato.cpp
393b17a3c8bbfeatures/netcdf-mpi

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"
/* -------------------------------------------------------------------------- */
__BEGIN_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();
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();
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();
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) {
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();
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;
}
__END_TAMAAS__
c4science · Help