Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/bem/bem_grid_polonski.cpp
84326222c0e9features/netcdf-mpi

Last Change
Raw File

bem_grid_polonski.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_polonski.hh"
#include "types.hh"
#include <limits>
#include <iomanip>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
BemGridPolonski::BemGridPolonski(Surface<Real> & surface):
BemGridKato(surface),
search_direction(surface.sizes(), 3),
projected_search_direction(surface.sizes(), 3)
{}
/* -------------------------------------------------------------------------- */
BemGridPolonski::~BemGridPolonski()
{}
/* -------------------------------------------------------------------------- */
Real BemGridPolonski::computeOptimalStep() {
this->applyInfluenceFunctions(this->search_direction, this->projected_search_direction);
centerOnContactArea(this->projected_search_direction);
Real num = 0;
Real den = 0;
const UInt N = search_direction.getNbPoints();
VectorProxy<Real>
p(nullptr, 3),
g(nullptr, 3),
r(nullptr, 3),
t(nullptr, 3);
UInt nc = search_direction.getNbComponents();
#pragma omp parallel for reduction(+:num, den) firstprivate(p, g, r, t)
for (UInt i = 0 ; i < N ; i++) {
p.setPointer(this->tractions.getInternalData() + i * nc);
g.setPointer(this->gaps.getInternalData() + i * nc);
r.setPointer(this->projected_search_direction.getInternalData() + i * nc);
t.setPointer(this->search_direction.getInternalData() + i * nc);
if (p(2) > 0) {
num += g.dot(t);
den += r.dot(t);
}
}
Real tau = num / den;
return tau;
}
/* -------------------------------------------------------------------------- */
void BemGridPolonski::updateSearchDirection(Real delta, Real G, Real G_old) {
Real * t_start = search_direction.getInternalData(),
* g_start = this->gaps.getInternalData(),
* p_start = this->tractions.getInternalData();
VectorProxy<Real>
t(nullptr, 3),
p(nullptr, 3),
g(nullptr, 3);
UInt nb_components = 3;
const UInt N = search_direction.getNbPoints();
#pragma omp parallel for firstprivate(t, p, g)
for (UInt i = 0 ; i < N ; i++) {
t.setPointer(t_start + i * nb_components);
g.setPointer(g_start + i * nb_components);
p.setPointer(p_start + i * nb_components);
if (p(2) > 0) {
t *= delta * G / G_old;
t += g;
} else {
t = 0;
}
}
}
/* -------------------------------------------------------------------------- */
Real BemGridPolonski::computeGradientNorm() {
Real G = 0;
VectorProxy<Real>
g(nullptr, 3),
p(nullptr, 3);
const UInt N = this->gaps.getNbPoints();
#pragma omp parallel for reduction(+:G) firstprivate(g, p)
for (UInt i = 0 ; i < N ; i++) {
g.setPointer(this->gaps.getInternalData() + i * g.dataSize());
p.setPointer(this->tractions.getInternalData() + i * p.dataSize());
if (p(2) > 0)
G += g.dot(g);
}
return G;
}
/* -------------------------------------------------------------------------- */
void BemGridPolonski::centerOnContactArea(Grid<Real, 2> &field) {
// Compute average values on contact
Real f0 = 0, f1 = 0, f2 = 0;
Real *f = field.getInternalData(), *p_start = this->tractions.getInternalData();
VectorProxy<Real> v(nullptr, 3);
VectorProxy<Real> p(nullptr, 3);
const UInt N = field.getNbPoints();
UInt nb_contact = 0;
#pragma omp parallel firstprivate(v, p)
{
#pragma omp for reduction(+:f0, f1, f2, nb_contact)
for (UInt i = 0 ; i < N ; i++) {
v.setPointer(f + i * v.dataSize());
p.setPointer(p_start + i * p.dataSize());
if (p(2) > 0) {
f0 += v(0);
f1 += v(1);
f2 += v(2);
nb_contact++;
}
}
#pragma omp single
{
f0 /= nb_contact; f1 /= nb_contact; f2 /= nb_contact;
}
#pragma omp barrier
#pragma omp for
for (UInt i = 0 ; i < N ; i++) {
v.setPointer(f + i * v.dataSize());
v(0) -= f0;
v(1) -= f1;
v(2) -= f2;
}
} // end #pragma omp parallel
}
/* -------------------------------------------------------------------------- */
Real BemGridPolonski::updateTractions(Real alpha) {
Real
*p_start = this->tractions.getInternalData(),
*g_start = this->gaps.getInternalData(),
*t_start = search_direction.getInternalData();
VectorProxy<Real>
p(nullptr, 3),
g(nullptr, 3),
t(nullptr, 3);
const UInt N = this->tractions.getNbPoints();
const UInt nb_components = this->tractions.getNbComponents();
UInt nb_iol = 0;
#pragma omp parallel firstprivate(p, g, t)
{
Real alpha_t_data[3] = {0};
VectorProxy<Real> alpha_t(alpha_t_data, 3);
#pragma omp for reduction(+:nb_iol)
for (UInt i = 0 ; i < N ; i++) {
p.setPointer(p_start + i * nb_components);
g.setPointer(g_start + i * nb_components);
t.setPointer(t_start + i * nb_components);
if (p(2) > 0) {
alpha_t = t;
alpha_t *= alpha;
p -= alpha_t;
BemGrid::projectOnFrictionCone(p, mu);
}
if (p(2) == 0 && lambda(i) < 0) {
nb_iol++;
p(2) -= alpha * lambda(i);
}
}
} // end #pragma omp parallel
return (Real) !(nb_iol > 0);
}
/* -------------------------------------------------------------------------- */
/// p_target is just a 3 component single value vector !
void BemGridPolonski::computeEquilibrium(Real tol, const Grid<Real, 1> & p_target) {
this->computeInfluence();
Real f = 0;
UInt n_iter = 0;
TAMAAS_ASSERT(p_target.dataSize() == 3, "Applied pressure must have only 3 components");
Real
delta = 0,
G = 0,
G_old = 1,
tau = 0;
setPressure(p_target);
// main loop
do {
this->computeDisplacementsFromTractions();
this->computeGaps(this->displacements);
centerOnContactArea(this->gaps);
beta = computeLambda(); // note: beta is member variable
G = computeGradientNorm();
// if (G == 0) break;
updateSearchDirection(delta, G, G_old);
G_old = G;
tau = computeOptimalStep();
delta = updateTractions(tau);
//Real f_ = 0;
//do {
// enforcePressureBalance(p_target);
// enforcePressureConstraints();
// f_ = computeEquilibriumError(p_target);
//} while (f_ > 1e-6);
lambda -= beta;
f = computeStoppingCriterion();
printFriendlyMessage(n_iter, f);
} while (f > tol && n_iter++ < this->max_iter);
}
/* -------------------------------------------------------------------------- */
void BemGridPolonski::setPressure(const Grid<Real, 1> & p_target) {
VectorProxy<Real> p(nullptr, 3);
const UInt N = this->tractions.getNbPoints();
Real *p_start = this->tractions.getInternalData();
#pragma omp parallel for firstprivate(p)
for (UInt i = 0 ; i < N ; i++) {
p.setPointer(p_start + i * p.dataSize());
p = p_target;
}
}
__END_TAMAAS__
c4science · Help