Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85468080
bem_gigipol.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, Sep 29, 11:42
Size
22 KB
Mime Type
text/x-c
Expires
Tue, Oct 1, 11:42 (2 d)
Engine
blob
Format
Raw Data
Handle
21164244
Attached To
rTAMAAS tamaas
bem_gigipol.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_gigipol.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
void BemGigipol::computeTractionsFromDisplacements() {
this->applyInverseInfluenceFunctions(this->true_displacements,this->surface_tractions);
}
Real BemGigipol::computeEquilibrium(Real epsilon,
Real mean_displacement) {
this->computeSpectralInfluenceOverDisplacement();
this->surface_t = 0.;
this->surface_r = 0.;
this->search_direction = 0.;
this->pold = 0.;
this->surface_displacements = 0.;
Real delta = 0.;
Real Gold = 1.;
Real f = 1e300;
Real fPrevious = 1e300;
UInt n = surface.size();
UInt size = n*n;
Real current_disp = SurfaceStatistics::computeAverage(true_displacements,0);
std::cout << "current disp "<<current_disp <<std::endl;
if (current_disp <= 0.) {
true_displacements = mean_displacement;
std::cout << "je re initialise "<< std::endl;
}
this->computeGaps();
this->optimizeToMeanDisplacement(mean_displacement);
std::ofstream file ("output.txt");
this->computeGaps();
convergence_iterations.clear();
nb_iterations = 0;max_iterations=10000;
while(f > epsilon && nb_iterations < max_iterations) {
fPrevious = f;
this->functional->computeGradFU();
this->search_direction= this->functional->getGradF();
Real gbar = this->computeMeanPressuresInNonContact();
this->search_direction -= gbar;
Real G = this->computeG();
this->updateT(G,Gold,delta);
Real tau = this->computeTau();
this->old_displacements=this->true_displacements;
Gold = G;
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
this->true_displacements(i)-=tau*this->surface_t(i);
}
//Projection on admissible space
this->computeGaps();
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(i)<0){
this->true_displacements(i)=this->surface(i);}
}
this->computeGaps();
delta = this->updateDisplacements(tau, mean_displacement);
//Projection on admissible space
this->computeGaps();
this->enforceMeanDisplacement(mean_displacement);
this->computeGaps();
this->computePressures();
f=this->computeStoppingCriterion();
if (nb_iterations % dump_freq ==0){
std::cout << "G vaut "<< G<< std::endl;
std::cout << "f vaut "<< f << std::endl;
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << " " << f-fPrevious << " "<< 'G' << " " << G << " "
<< std::endl;
UInt n = surface.size();
UInt size = n*n;
Real crit=0;
Real disp_norm=0;
for (UInt i = 0; i < size; ++i) {
crit += this->search_direction(i)*this->search_direction(i);
disp_norm += (true_displacements(i)*true_displacements(i));
}
file<< std::scientific << std::setprecision(10)<< nb_iterations << " "
<<crit/disp_norm<< " "<<f<< std::endl;
}
convergence_iterations.push_back(f);
++nb_iterations;
}
this->computePressures();
return f;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeEquilibrium2(Real epsilon,
Real mean_pressure) {
this->computeSpectralInfluenceOverDisplacement();
this->surface_t = 0.;
this->surface_r = 0.;
this->search_direction = 0.;
this->pold = 0.;
this->surface_displacements = 0.;
Real delta = 0.;
Real Gold = 1.;
Real f = 1e300;
Real fPrevious = 1e300;
UInt n = surface.size();
UInt size = n*n;
Real current_disp = SurfaceStatistics::computeAverage(surface,0.);
true_displacements = current_disp;
this->computeGaps();
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(i)<0){
this->true_displacements(i)=this->surface(i);}
}
this->computeGaps();
// this->optimizeToMeanDisplacement(mean_displacement);
std::ofstream file ("output.txt");
convergence_iterations.clear();
nb_iterations = 0;max_iterations=200000;
while(f > epsilon && nb_iterations < max_iterations) {
fPrevious = f;
this->functional->computeGradFU();
this->search_direction= this->functional->getGradF();
Real gbar = this->computeMeanPressuresInNonContact();
this->search_direction += (2*mean_pressure+gbar);
Real G = this->computeG();
this->updateT(G,Gold,delta);
Real tau = this->computeTau(mean_pressure);
this->old_displacements=this->true_displacements;
Gold = G;
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
this->true_displacements(i)-=tau*this->surface_t(i);
}
//Projection on admissible space
this->computeGaps();
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(i)<0){
this->true_displacements(i)=this->surface(i);}
}
this->computeGaps();
delta = this->updateDisplacements(tau, 0.);
//Projection on admissible space
this->computeGaps();
this->computePressures();
f=this->computeStoppingCriterion();
this->computeTractionsFromDisplacements();
this->functional->computeGradFU();
const Surface<Real> & gradF = this->functional->getGradF();
Real min = SurfaceStatistics::computeMinimum(gradF);
if (f< epsilon && epsilon<std::abs(min+mean_pressure)) f=3.;
if (nb_iterations % dump_freq ==0){
std::cout << "G vaut "<< G<< std::endl;
std::cout << "f vaut "<< f << std::endl;
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << " " << f-fPrevious << " "<< 'G' << " " << G << " "
<< std::endl;
UInt n = surface.size();
UInt size = n*n;
Real crit=0;
Real disp_norm=0;
for (UInt i = 0; i < size; ++i) {
crit += this->search_direction(i)*this->search_direction(i);
disp_norm += (true_displacements(i)*true_displacements(i));
}
file<< std::scientific << std::setprecision(10)<< nb_iterations << " "
<<crit/disp_norm<< " "<<f<< std::endl;
}
convergence_iterations.push_back(f);
++nb_iterations;
}
this->computePressures();
return f;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeEquilibrium2init(Real epsilon,
Real mean_pressure,Surface<Real> & init) {
this->computeSpectralInfluenceOverDisplacement();
this->surface_t = 0.;
this->surface_r = 0.;
this->search_direction = 0.;
this->pold = 0.;
this->surface_displacements = 0.;
Real delta = 0.;
Real Gold = 1.;
Real f = 1e300;
Real fPrevious = 1e300;
UInt n = surface.size();
UInt size = n*n;
this->true_displacements=init;
this->computeGaps();
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(i)<0){
this->true_displacements(i)=this->surface(i);}
}
this->computeGaps();
std::ofstream file ("output.txt");
convergence_iterations.clear();
nb_iterations = 0;max_iterations=200000;
while(f > epsilon && nb_iterations < max_iterations) {
fPrevious = f;
this->functional->computeGradFU();
this->search_direction= this->functional->getGradF();
Real gbar = this->computeMeanPressuresInNonContact();
this->search_direction += (2*mean_pressure+gbar);
Real G = this->computeG();
this->updateT(G,Gold,delta);
Real tau = this->computeTau(mean_pressure);
tau=0.01*tau;
this->old_displacements=this->true_displacements;
Gold = G;
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
this->true_displacements(i)-=tau*this->surface_t(i);
}
//Projection on admissible space
this->computeGaps();
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(i)<0){
this->true_displacements(i)=this->surface(i);}
}
this->computeGaps();
delta = this->updateDisplacements(tau, 0.);
//Projection on admissible space
this->computeGaps();
this->computePressures();
f=this->computeStoppingCriterion();
this->computeTractionsFromDisplacements();
this->functional->computeGradFU();
const Surface<Real> & gradF = this->functional->getGradF();
Real min = SurfaceStatistics::computeMinimum(gradF);
if (f< epsilon && epsilon<std::abs(min+mean_pressure)) f=3.;
if (nb_iterations % dump_freq ==0){
std::cout << "G vaut "<< G<< std::endl;
std::cout << "f vaut "<< f << std::endl;
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << " " << f-fPrevious << " "<< 'G' << " " << G << " "
<< std::endl;
UInt n = surface.size();
UInt size = n*n;
Real crit=0;
Real disp_norm=0;
for (UInt i = 0; i < size; ++i) {
crit += this->search_direction(i)*this->search_direction(i);
disp_norm += (true_displacements(i)*true_displacements(i));
}
file<< std::scientific << std::setprecision(10)<< nb_iterations << " "
<<crit/disp_norm<< " "<<f<< std::endl;
}
convergence_iterations.push_back(f);
++nb_iterations;
}
this->computePressures();
return f;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeEquilibriuminit(Real epsilon,
Real mean_displacement,Surface<Real> & init) {
this->computeSpectralInfluenceOverDisplacement();
this->surface_t = 0.;
this->surface_r = 0.;
this->search_direction = 0.;
this->pold = 0.;
this->surface_displacements = 0.;
Real delta = 0.;
Real Gold = 1.;
Real f = 1e300;
Real fPrevious = 1e300;
UInt n = surface.size();
UInt size = n*n;
this->true_displacements=init;
Real current_disp = SurfaceStatistics::computeAverage( this->true_displacements,0);
std::cout << "current disp "<<current_disp <<std::endl;
if (current_disp == 0.) {
true_displacements = mean_displacement;
std::cout << "je reinitialise "<<current_disp <<std::endl;
}
this->computeGaps();
this->optimizeToMeanDisplacement(mean_displacement);
std::ofstream file ("output.txt");
this->computeGaps();
convergence_iterations.clear();
nb_iterations = 0;max_iterations=100000;
while(f > epsilon && nb_iterations < max_iterations) {
fPrevious = f;
this->functional->computeGradFU();
this->search_direction= this->functional->getGradF();
Real gbar = this->computeMeanPressuresInNonContact();
this->search_direction -= gbar;
Real G = this->computeG();
this->updateT(G,Gold,delta);
Real tau = this->computeTau();
this->old_displacements=this->true_displacements;
Gold = G;
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
this->true_displacements(i)-=tau*this->surface_t(i);
}
//Projection on admissible space
this->computeGaps();
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(i)<0){
this->true_displacements(i)=this->surface(i);}
}
this->computeGaps();
delta = this->updateDisplacements(tau, mean_displacement);
//Projection on admissible space
this->computeGaps();
this->enforceMeanDisplacement(mean_displacement);
this->computeGaps();
this->computePressures();
f=this->computeStoppingCriterion();
if (nb_iterations % dump_freq ==0){
std::cout << "G vaut "<< G<< std::endl;
std::cout << "f vaut "<< f << std::endl;
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << " " << f-fPrevious << " "<< 'G' << " " << G << " "
<< std::endl;
UInt n = surface.size();
UInt size = n*n;
Real crit=0;
Real disp_norm=0;
for (UInt i = 0; i < size; ++i) {
crit += this->search_direction(i)*this->search_direction(i);
disp_norm += (true_displacements(i)*true_displacements(i));
}
file<< std::scientific << std::setprecision(10)<< nb_iterations << " "
<<crit/disp_norm<< " "<<f<< std::endl;
}
convergence_iterations.push_back(f);
++nb_iterations;
}
this->computePressures();
return f;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeStoppingCriterion() {
UInt n = surface.size();
UInt size = n*n;
Real res = 0;
Real t_sum = std::abs(SurfaceStatistics::computeSum(this->true_displacements));
this->computeTractionsFromDisplacements();
this->functional->computeGradFU();
const Surface<Real> & gradF = this->functional->getGradF();
Real min = SurfaceStatistics::computeMinimum(gradF);
#pragma omp parallel for reduction(+:res)
for (UInt i = 0; i < size; ++i) {
res += std::abs((gradF(i)-min)* ( 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);
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::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->gap(i) > 0.) {
Real val = this->search_direction(i);
res += val*val;}
}
STOPTIMER("computeG");
return res;
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeMeanPressuresInNonContact(){
STARTTIMER("computeMeanPressuresInNonContact");
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->gap(i) == 0. ) continue;
++nb_contact;
res += this->search_direction(i);
}
res /= nb_contact;
STOPTIMER("computeMeanPressuresInNonContact");
return res;
}
/* -------------------------------------------------------------------------- */
void BemGigipol::optimizeToMeanDisplacement(Real imposed_mean) {
Real target_value = imposed_mean - SurfaceStatistics::computeAverage(surface, 0);
UInt n = surface.size();
UInt size = n * n;
// Initial guesses for upper and lower bound
Real step_min = -10;
Real step_max = 10;
// Gaps for upper and lower bound
Real gap_min = positiveGapAverage(step_min);
Real gap_max = positiveGapAverage(step_max);
UInt max_expansion = 8;
// Expand bounds if necessary
for (UInt i = 0 ; gap_max < target_value && i < max_expansion ; i++, step_max *= 10)
gap_max = positiveGapAverage(step_max);
for (UInt i = 0 ; gap_min > target_value && i < max_expansion ; i++, step_min *= 10)
gap_min = positiveGapAverage(step_min);
Real g = 0.;
Real epsilon = 1e-12;
Real step = 0;
while (fabs(g - target_value) > epsilon) {
step = (step_min + step_max) / 2.;
g = positiveGapAverage(step);
if (g > target_value) step_max = step;
else if (g < target_value) step_min = step;
else {
step_max = step;
step_min = step;
}
}
step = (step_min + step_max) / 2.;
#pragma omp parallel for
for (UInt i = 0 ; i < size ; i++) {
gap(i) += step;
if (gap(i) < 0) gap(i) = 0;
true_displacements(i) = gap(i) + surface(i);
}
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::positiveGapAverage(Real shift) {
UInt n = surface.size();
Real res = 0;
#pragma omp parallel for reduction(+: res)
for (UInt i = 0 ; i < n * n ; i++) {
Real shifted_gap = gap(i) + shift;
res += shifted_gap * (shifted_gap > 0);
}
return res / (n * n);
}
/* -------------------------------------------------------------------------- */
void BemGigipol::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->gap(i) == 0.) this->surface_t(i) = 0.;
else{
this->surface_t(i) *= factor;
this->surface_t(i) += this->search_direction(i);
}
}
STOPTIMER("updateT");
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeTau() {
STARTTIMER("computeOptimalStep");
const UInt n = surface.size();
const UInt size = n*n;
this->applyInverseInfluenceFunctions(surface_t, surface_r);
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->gap(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->gap(i) == 0) continue;
tau_sum1 += this->search_direction(i)*surface_t(i);
tau_sum2 += surface_r(i)*surface_t(i);
}
Real tau = tau_sum1/tau_sum2;
STOPTIMER("computeTau");
return tau;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeTau(Real mean_pressure) {
STARTTIMER("computeOptimalStep");
const UInt n = surface.size();
const UInt size = n*n;
this->applyInverseInfluenceFunctions(surface_t, surface_r);
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->gap(i) == 0) continue;
++nb_contact;
rbar += surface_r(i);
}
rbar /= nb_contact;
surface_r += 2*mean_pressure+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->gap(i) == 0) continue;
tau_sum1 += this->search_direction(i)*surface_t(i);
tau_sum2 += surface_r(i)*surface_t(i);
}
Real tau = tau_sum1/tau_sum2;
STOPTIMER("computeTau");
return tau;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::updateDisplacements(Real tau, Real mean_displacement) {
STARTTIMER("updateDisplacements");
unsigned int n = surface.size();
unsigned int size = n*n;
//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->search_direction(i) <0 && this->gap(i) == 0.){
this->true_displacements(i) -= tau* this->search_direction(i);
++nb_iol;
}
}
Real delta = 0;
if (nb_iol > 0) delta = 0.;
else delta = 1.;
return delta;
STOPTIMER("updateDisplacements");
}
/* -------------------------------------------------------------------------- */
void BemGigipol::enforceMeanDisplacement(Real mean_displacement) {
STARTTIMER("enforceMeanDisplacement");
UInt n = surface.size();
UInt size = n*n;
Real moyenne_surface=SurfaceStatistics::computeAverage(this->surface, 0);
Real average = SurfaceStatistics::computeAverage(this->true_displacements, 0);
Real factor = (mean_displacement-moyenne_surface) / (average-moyenne_surface);
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->true_displacements(i) = factor*(this->true_displacements(i)-this->surface(i)) +this->surface(i);
}
STOPTIMER("enforceMeanDisplacement");
}
/* -------------------------------------------------------------------------- */
void BemGigipol::computePressures() {
this->computeTractionsFromDisplacements();
this->functional->computeGradFU();
const Surface<Real> & gradF = this->functional->getGradF();
Real min = SurfaceStatistics::computeMinimum(gradF);
this->surface_tractions -= min;
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
Event Timeline
Log In to Comment