Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/bem/bem_kato.cpp
fdea6b4ef4a4features/declarative

Last Change
Raw File

bem_kato.cpp
View Options

/**
*
* @author Guillaume Anciaux <guillaume.anciaux@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_kato.hh"
#include "surface.hh"
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <vector>
/* -------------------------------------------------------------------------- */
namespace tamaas {
Real BemKato::linescan(Real scale, Real pressure) {
updatePressure(scale);
shiftPressure(pressure);
truncatePressure();
Real res = computeF();
return res;
}
/* -------------------------------------------------------------------------- */
Real BemKato::linesearch(Real& hmax, Real fold, Real pressure,
int search_flag) {
if (search_flag) {
Real h = hmax;
// Real fold = bem.computeF();
// if (fold == 0) fold = 1e300;
Real f = linescan(h, pressure);
if (f < fold)
return f;
while (f > fold) {
h *= 0.5;
if (h < 1e-3)
throw 1;
f = linescan(h, pressure);
}
f = linescan(h, pressure);
// if (hmax / h > 10) hmax /=2;
return f;
}
return linescan(hmax, pressure);
}
/* -------------------------------------------------------------------------- */
Real BemKato::computeEquilibrium(Real epsilon, Real pressure) {
// UInt n = surface.size();
// UInt size = n*n;
this->computeSpectralInfluenceOverDisplacement();
this->computeDisplacementsFromTractions();
this->functional->computeGradFP();
this->backupTractions();
Real f = 1.;
Real fPrevious = 1e300;
this->nb_iterations = 0;
this->convergence_iterations.clear();
while (f > epsilon && this->nb_iterations < this->max_iterations) {
fPrevious = f;
this->computeDisplacementsFromTractions();
this->functional->computeGradFP();
this->backupTractions();
try {
f = linescan(1., pressure);
} catch (int e) {
std::cout << " line search failed " << std::endl;
f = linescan(1, pressure);
nb_iterations = 0;
break;
}
if (nb_iterations % dump_freq == 0) {
// std::cout << std::scientific << std::setprecision(10) <<
// nb_iterations << " " << f << " " << f-fold << " " <<
// ((f-fold)/forigin) << std::endl;
std::cout << std::scientific << std::setprecision(10) << nb_iterations
<< " " << f << " " << f - fPrevious << std::endl;
}
convergence_iterations.push_back(f);
++nb_iterations;
}
this->computeTrueDisplacements();
this->computeGaps();
return f;
}
/* -------------------------------------------------------------------------- */
void BemKato::updatePressure(Real scale) {
UInt n = surface.size();
UInt size = n * n;
const Surface<Real>& gradF = this->functional->getGradF();
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->surface_tractions(i) =
this->surface_tractions_backup(i) - gradF(i) * scale;
}
}
/* -------------------------------------------------------------------------- */
void BemKato::backupTractions() {
UInt n = surface.size();
UInt size = n * n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->surface_tractions_backup(i) = this->surface_tractions(i);
}
}
/* -------------------------------------------------------------------------- */
Real BemKato::positivePressure(Real step) {
UInt n = surface.size();
UInt size = n * n;
Real p_tot = 0.0;
#pragma omp parallel for reduction(+ : p_tot)
for (UInt i = 0; i < size; ++i) {
Real sh_press = this->surface_tractions(i) + step;
if (sh_press > max_pressure)
p_tot += max_pressure;
else
p_tot += sh_press * (sh_press > 0);
}
return p_tot / n / n;
}
/* -------------------------------------------------------------------------- */
void BemKato::shiftPressure(Real required_pressure) {
Real step_min = -10;
Real step_max = 10;
Real p_max = positivePressure(step_max);
Real p_min = positivePressure(step_min);
for (UInt i = 0; p_max < required_pressure && i < 8; ++i, step_max *= 10) {
p_max = positivePressure(step_max);
}
for (UInt i = 0; p_min > required_pressure && i < 8; ++i, step_min *= 10) {
p_min = positivePressure(step_min);
}
Real p = positivePressure(0.0);
Real epsilon = 1e-12;
while (std::abs(step_min - step_max) > epsilon) {
Real step = (step_min + step_max) / 2;
p = positivePressure(step);
if (p > required_pressure)
step_max = step;
else if (p < required_pressure)
step_min = step;
else {
step_max = step;
step_min = step;
}
}
UInt n = surface.size();
UInt size = n * n;
// shift the pressure so that satisfies the constraint
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->surface_tractions(i) += (step_max + step_min) / 2;
}
}
/* -------------------------------------------------------------------------- */
void BemKato::truncatePressure() {
UInt n = surface.size();
UInt size = n * n;
// shift the pressure so that satisfies the constraint
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) > max_pressure)
this->surface_tractions(i) = max_pressure;
else
this->surface_tractions(i) *= (this->surface_tractions(i) > 0);
}
}
/* -------------------------------------------------------------------------- */
void BemKato::computeTrueDisplacements() {
UInt n = surface.size();
UInt size = n * n;
Real shift = 1e300;
#pragma omp parallel for reduction(min : shift)
for (UInt i = 0; i < size; i++) {
if (surface_displacements(i) - shift - surface(i) < 0. &&
surface_tractions(i) < max_pressure) {
shift = surface_displacements(i) - surface(i);
}
}
this->true_displacements = surface_displacements;
this->true_displacements -= shift;
}
/* -------------------------------------------------------------------------- */
Real BemKato::computeF() {
UInt n = surface.size();
UInt size = n * n;
Real res = 0;
Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions));
computeTrueDisplacements();
#pragma omp parallel for reduction(+ : res)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) == this->max_pressure)
continue;
res += std::abs(surface_tractions(i) *
(this->true_displacements(i) - surface(i)));
// res +=
// this->surface_tractions[i].real()
// *(surface_displacements[i].real() - surface[i].real());
}
return res / std::abs(t_sum);
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas
c4science · Help