Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F69243240
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
Sun, Jun 30, 22:05
Size
8 KB
Mime Type
text/x-c
Expires
Tue, Jul 2, 22:05 (2 d)
Engine
blob
Format
Raw Data
Handle
18686990
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 <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeEquilibrium(Real epsilon,
Real pressure){
this->computeSpectralInfluenceOverDisplacement();
this->surface_t = 0.;
this->surface_r = 0.;
this->surface_spectral_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->computeDisplacements();
this->computeGaps();
Real gbar = this->computeMeanGapsInContact();
this->gap -= gbar;
Real G = this->computeG();
std::cout << "G vaut "<< G<< std::endl;
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);
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::endl;
}
convergence_iterations.push_back(f);
++nb_iterations;
}
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");
unsigned int n = surface.size();
unsigned int size = n*n;
Real res = 0.;
UInt nb_contact = 0;
#pragma omp parallel for reduction(+: nb_contact,res)
for (unsigned int i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) continue;
++nb_contact;
res += this->gap(i);
}
res /= nb_contact;
STOPTIMER("computeMeanGapsInContact");
return res;
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeG(){
STARTTIMER("computeG");
unsigned int n = surface.size();
unsigned int size = n*n;
Real res = 0.;
#pragma omp parallel for reduction(+: res)
for (unsigned int i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) continue;
Real val = this->gap(i);
res += val*val;
}
STOPTIMER("computeG");
return res;
}
/* -------------------------------------------------------------------------- */
void BemPolonski::updateT(Real G, Real Gold, Real delta){
STARTTIMER("updateT");
unsigned int n = surface.size();
unsigned int size = n*n;
Real factor = delta*G/Gold;
#pragma omp parallel for
for (unsigned int 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");
surface_spectral_r = surface_t;
surface_spectral_r.FFTTransform(nthreads);
unsigned int n = surface.size();
unsigned int size = n*n;
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
surface_spectral_r(i) *= this->surface_spectral_influence_disp(i);
}
surface_spectral_r.FFTITransform(surface_r,nthreads);
Real rbar = 0;
UInt nb_contact = 0;
#pragma omp parallel for reduction(+: nb_contact,rbar)
for (unsigned int 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 (unsigned int 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");
unsigned int n = surface.size();
unsigned int size = n*n;
#pragma omp parallel for
for (unsigned int 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 (unsigned int 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;
}
}
Real delta = 0;
if (nb_iol > 0) delta = 0.;
else delta = 1.;
STOPTIMER("updateTractions");
return delta;
}
/* -------------------------------------------------------------------------- */
void BemPolonski::enforcePressureBalance(Real applied_pressure){
STARTTIMER("enforcePressureBalance");
unsigned int n = surface.size();
unsigned int size = n*n;
Real pressure = 0.;
#pragma omp parallel for reduction(+: pressure)
for (unsigned int i = 0; i < size; ++i) {
pressure += this->surface_tractions(i);
}
pressure *= 1./size;
for (unsigned int i = 0; i < size; ++i) {
this->surface_tractions(i) *= applied_pressure/pressure;
}
STOPTIMER("enforcePressureBalance");
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeF() {
UInt n = surface.size();
UInt size = n*n;
Real res = 0;
Real d_max = SurfaceStatistics::computeMaximum(surface_displacements);
Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions));
Real s_max = SurfaceStatistics::computeMaximum(surface);
#pragma omp parallel for reduction(+:res)
for (UInt i = 0; i < size; ++i) {
res += std::abs(surface_tractions(i) * (surface_displacements(i)-d_max-surface(i)+s_max));
// res +=
// this->surface_tractions[i].real()
// *(surface_displacements[i].real() - surface[i].real());
}
return res/std::abs(t_sum);
}
/* -------------------------------------------------------------------------- */
Event Timeline
Log In to Comment