Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/bem_gigipol.cpp
f05be88c05d7maxwell_viscoelastic…

Last Change
Raw File

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__
c4science · Help