Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65894503
bem_polonski.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Jun 6, 22:32
Size
15 KB
Mime Type
text/x-c
Expires
Sat, Jun 8, 22:32 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18140524
Attached To
rTAMAAS tamaas
bem_polonski.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_polonski.hh"
#include "fft_plan_manager.hh"
#include "fftransform.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
BemPolonski::BemPolonski(Surface<Real> & p) :
BemFFTBase(p),
surface_t(p.size(),p.getL()),
surface_r(p.size(),p.getL()),
pold(p.size(),p.getL()),
p_sat(-1.){
e_star = 1.;
max_iterations = 2000;
}
/* -------------------------------------------------------------------------- */
BemPolonski::~BemPolonski() {}
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeEquilibrium(Real epsilon,
Real pressure) {
this->computeSpectralInfluenceOverDisplacement();
this->surface_t = 0.;
this->surface_r = 0.;
this->pold = 0.;
this->surface_displacements = 0.;
Real delta = 0.;
Real Gold = 1.;
Real f = 1e300;
Real fPrevious = 1e300;
Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions);
if (current_pressure <= 0.) surface_tractions = pressure;
this->enforcePressureBalance(pressure);
convergence_iterations.clear();
nb_iterations = 0;
while(f > epsilon && nb_iterations < max_iterations) {
fPrevious = f;
this->computeDisplacementsFromTractions();
this->computeGaps();
Real gbar = this->computeMeanGapsInContact();
this->gap -= gbar;
Real G = this->computeG();
this->updateT(G,Gold,delta);
Real tau = this->computeTau();
pold = this->surface_tractions;
Gold = G;
delta = this->updateTractions(tau);
//Projection on admissible space
this->enforcePressureBalance(pressure);
// std::cout << std::scientific << std::setprecision(10)
// << SurfaceStatistics::computeAverage(surface_tractions) << std::fixed << std::endl;
//
f = this->computeF();
if (nb_iterations % dump_freq == 0){
Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions);
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << " " << f-fPrevious << " " << A
<< std::fixed << std::endl;
}
convergence_iterations.push_back(f);
++nb_iterations;
}
this->computeTrueDisplacements();
return f;
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeEquilibriuminit(Real epsilon,
Real pressure, Surface<Real> & init) {
this->computeSpectralInfluenceOverDisplacement();
this->surface_t = 0.;
this->surface_r = 0.;
this->pold = 0.;
this->surface_displacements = 0.;
Real delta = 0.;
Real Gold = 1.;
Real f = 1e300;
Real fPrevious = 1e300;
this->surface_tractions=init;
Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions);
if (current_pressure <= 0.) surface_tractions = pressure;
this->enforcePressureBalance(pressure);
convergence_iterations.clear();
nb_iterations = 0;
this->computeDisplacementsFromTractions();
f=this->computeF();
std::cout << std::scientific << std::setprecision(10)
<<" "
<< f << std::endl;
while(f > epsilon && nb_iterations < max_iterations) {
fPrevious = f;
this->computeDisplacementsFromTractions();
this->computeGaps();
Real gbar = this->computeMeanGapsInContact();
this->gap -= gbar;
Real G = this->computeG();
this->updateT(G,Gold,delta);
Real tau = this->computeTau();
pold = this->surface_tractions;
Gold = G;
delta = this->updateTractions(tau);
//Projection on admissible space
this->enforcePressureBalance(pressure);
// std::cout << std::scientific << std::setprecision(10)
// << SurfaceStatistics::computeAverage(surface_tractions) << std::fixed << std::endl;
//
f = this->computeF();
if (nb_iterations % dump_freq == 0){
Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions);
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << " " << f-fPrevious << " " << A
<< std::fixed << std::endl;
}
convergence_iterations.push_back(f);
++nb_iterations;
}
this->computeTrueDisplacements();
return f;
}
/* -------------------------------------------------------------------------- */
void BemPolonski::computeGaps() {
UInt size = surface.size();
#pragma omp parallel for
for (UInt i=0; i < size*size; i++) {
gap(i) = surface_displacements(i) - surface(i);
}
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeMeanGapsInContact(){
STARTTIMER("computeMeanGapsInContact");
UInt n = surface.size();
UInt size = n*n;
Real res = 0.;
UInt nb_contact = 0;
#pragma omp parallel for reduction(+: nb_contact,res)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0 ||std::abs(this->surface_tractions(i) - p_sat)< 1.0e-10 ) continue;
++nb_contact;
res += this->gap(i);
}
res /= nb_contact;
STOPTIMER("computeMeanGapsInContact");
return res;
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeG(){
STARTTIMER("computeG");
UInt n = surface.size();
UInt size = n*n;
Real res = 0.;
#pragma omp parallel for reduction(+: res)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0 ||std::abs(this->surface_tractions(i) - p_sat)< 1.0e-10 ) continue;
Real val = this->gap(i);
res += val*val;
}
STOPTIMER("computeG");
return res;
}
/* -------------------------------------------------------------------------- */
void BemPolonski::updateT(Real G, Real Gold, Real delta){
STARTTIMER("updateT");
UInt n = surface.size();
UInt size = n*n;
Real factor = delta*G/Gold;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) this->surface_t(i) = 0.;
else{
this->surface_t(i) *= factor;
this->surface_t(i) += this->gap(i);
}
}
STOPTIMER("updateT");
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeTau(){
STARTTIMER("computeTau");
this->applyInfluenceFunctions(surface_t, surface_r);
const UInt size = surface_t.size() * surface_t.size();
Real rbar = 0;
UInt nb_contact = 0;
#pragma omp parallel for reduction(+: nb_contact,rbar)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) continue;
++nb_contact;
rbar += surface_r(i);
}
rbar /= nb_contact;
surface_r -= rbar;
Real tau_sum1 = 0.,tau_sum2 = 0.;
#pragma omp parallel for reduction(+: tau_sum1, tau_sum2)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) continue;
tau_sum1 += this->gap(i)*surface_t(i);
tau_sum2 += surface_r(i)*surface_t(i);
}
Real tau = tau_sum1/tau_sum2;
STOPTIMER("computeTau");
return tau;
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::updateTractions(Real tau){
STARTTIMER("updateTractions");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) continue;
this->surface_tractions(i) -= tau*this->surface_t(i);
if (this->surface_tractions(i) < 0.) this->surface_tractions(i) = 0;
}
//compute number of interpenetration without contact
UInt nb_iol = 0;
#pragma omp parallel for reduction(+: nb_iol)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0. &&
this->gap(i) < 0.) {
this->surface_tractions(i) -= tau*this->gap(i);
++nb_iol;
}
}
//compute number of interpenetration without contact
UInt nb_sl = 0;
#pragma omp parallel for reduction(+: nb_sl)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) ==p_sat &&
0.< this->gap(i) ) {
this->surface_tractions(i) -= tau*this->gap(i);
++nb_sl;
}
}
Real delta = 0;
if (nb_iol > 0 || nb_sl> 0) delta = 0.;
else delta = 1.;
STOPTIMER("updateTractions");
return delta;
}
/* -------------------------------------------------------------------------- */
void BemPolonski::enforcePressureBalance(Real applied_pressure){
STARTTIMER("enforcePressureBalance");
UInt n = surface.size();
UInt size = n*n;
// If we have no saturation, scale to match applied_pressure
if (this->p_sat <= 0. ||
SurfaceStatistics::computeMaximum(this->surface_tractions) < this->p_sat) {
Real pressure = 0.;
#pragma omp parallel for reduction(+: pressure)
for (UInt i = 0; i < size; ++i) {
pressure += this->surface_tractions(i);
}
pressure *= 1./size;
this->surface_tractions *= applied_pressure/pressure;
}
// If we have saturation, use secant method to find zero of saturationFunction
else {
// // Checking if we can find a zero
// Real limit =
// this->p_sat * SurfaceStatistics::computeContactArea(this->surface_tractions)
// - applied_pressure;
// if (limit < 0) SURFACE_FATAL("Saturation pressure is too small for applied pressure");
// // Initial points
// Real x_n_2 = 0., x_n_1 = 1., x_n = 0.;
// Real f_n_2 = 0., f_n_1 = 0., f_n = 0.;
// // Secant loop
// do {
// f_n_2 = saturationFunction(x_n_2, applied_pressure);
// f_n_1 = saturationFunction(x_n_1, applied_pressure);
// x_n = x_n_1 - f_n_1 * (x_n_1 - x_n_2) / (f_n_1 - f_n_2);
// f_n = saturationFunction(x_n, applied_pressure);
// x_n_2 = x_n_1;
// x_n_1 = x_n;
// } while (std::abs(f_n) > 1e-16);
// this->surface_tractions *= x_n;
// // Truncating pressures
//#pragma omp parallel for
// for (UInt i = 0 ; i < size ; i++) {
// if (this->surface_tractions(i) > this->p_sat)
// this->surface_tractions(i) = this->p_sat;
// }
Real mu_I_sat=0.;
#pragma omp parallel for reduction(+: mu_I_sat)
for (UInt i = 0; i < size; ++i) {
if ( this->p_sat< this->surface_tractions(i) ){
mu_I_sat+=1.;
}
}
Real pc=0.;
#pragma omp parallel for reduction(+: pc)
for (UInt i = 0; i < size; ++i) {
if ( this->surface_tractions(i) < this->p_sat ){
pc+=this->surface_tractions(i);}
}
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if ( this->surface_tractions(i) < this->p_sat ){
this->surface_tractions(i)= this->surface_tractions(i)*((size*applied_pressure-this->p_sat*mu_I_sat)/pc);}
}
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if ( this->p_sat< this->surface_tractions(i) ){
this->surface_tractions(i)=this->p_sat;
}
}
}
STOPTIMER("enforcePressureBalance");
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::saturationFunction(Real alpha, Real applied_pressure) {
const UInt n = surface.size();
const UInt size = n*n;
Real sum = 0.;
#pragma omp parallel for reduction(+: sum)
for (UInt i = 0 ; i < size ; i++) {
Real alpha_p = alpha * this->surface_tractions(i);
Real alpha_p_I = (alpha_p > this->p_sat)? alpha_p - this->p_sat : 0.;
sum += alpha_p - alpha_p_I;
}
sum /= size;
return sum - applied_pressure;
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::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->p_sat) 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);
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
void BemPolonski::computeTrueDisplacements(){
this->applyInfluenceFunctions(this->surface_tractions, this->surface_displacements);
this->true_displacements = this->surface_displacements;
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. &&( this->p_sat < 0 || this->surface_tractions(i) < this->p_sat )){
shift = surface_displacements(i) - surface(i);
}
}
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
true_displacements(i) = surface_displacements(i) - shift;
}
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
Event Timeline
Log In to Comment