Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/bem_kato.cpp
562a3d65fd33features/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 <vector>
#include "surface.hh"
#include "bem_kato.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_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){
STARTTIMER("updatePressure");
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;
}
STOPTIMER("updatePressure");
}
/* -------------------------------------------------------------------------- */
void BemKato::backupTractions(){
STARTTIMER("switchPressure");
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);
}
STOPTIMER("switchPressure");
}
/* -------------------------------------------------------------------------- */
Real BemKato::positivePressure(Real step){
STARTTIMER("positivePressure");
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);
}
STOPTIMER("positivePressure");
return p_tot/n/n;
}
/* -------------------------------------------------------------------------- */
void BemKato::shiftPressure(Real required_pressure){
Real step_min = -10;
Real step_max = 10;
STARTTIMER("shiftPressureInitialGuess");
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;
STOPTIMER("shiftPressureInitialGuess");
STARTTIMER("shiftPressureDichotomy");
while (fabs(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;
}
}
STOPTIMER("shiftPressureDichotomy");
STARTTIMER("shiftPressureTrueShift");
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;
}
STOPTIMER("shiftPressureTrueShift");
}
/* -------------------------------------------------------------------------- */
void BemKato::truncatePressure(){
STARTTIMER("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);
}
STOPTIMER("truncatePressure");
}
/* -------------------------------------------------------------------------- */
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);
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
c4science · Help