Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62399348
bem_gigi.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, May 12, 22:33
Size
11 KB
Mime Type
text/x-c
Expires
Tue, May 14, 22:33 (2 d)
Engine
blob
Format
Raw Data
Handle
17642754
Attached To
rTAMAAS tamaas
bem_gigi.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_gigi.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
Real BemGigi::computeEquilibrium(Real epsilon,
Real mean_displacement) {
this->computeSpectralInfluenceOverDisplacement();
Real Rold = 1.;
//Real Gold=1.;
Real f = 1e300;
this->search_direction = 0.;
this->true_displacements =0;
this->true_displacements = mean_displacement;
this->computeGaps();
this->optimizeToMeanDisplacement(mean_displacement);
this->computeGaps();
std::cout << "moyenne deplacement "<< SurfaceStatistics::computeAverage(this->true_displacements, 0)<< std::endl;
convergence_iterations.clear();
nb_iterations = 0;
UInt crit2=10;
max_iterations =1000;
while (f> epsilon && nb_iterations++ < max_iterations) {
this->computeGaps();
this->functional->computeGradF();
Real R = this->functional->computeGradientNorm();
std::cout << "norme gradient "<< R << std::endl;
this->updateSearchDirection(R / Rold);
Rold = R;
Real alpha = this->computeOptimalStep();
std::cout << "alpha vaut "<< alpha<< std::endl;
this->old_displacements=this->true_displacements;
this->updateDisplacements(alpha);
std::cout << "moyenne deplacement "<< SurfaceStatistics::computeAverage(this->true_displacements, 0)<< std::endl;
this->computeGaps();
this->optimizeToMeanDisplacement(mean_displacement); //espace admissible
this->computeGaps();
this->computePressures(mean_displacement); //p=K-1*u
Real G=this->computeG();
G=0.;
// std::cout << "G vaut "<< G<< std::endl;
UInt inner_loop_cpt=0;
//Inner loop
while(G > 0. && inner_loop_cpt < crit2 ) {
++inner_loop_cpt;
this->computeGaps();
this->functional->computeGradF();
this->updateSearchDirection(R / Rold);
Real tau = this->computeTau();
//Gold = G;
this->updateDisplacementsB(tau);
this->computeGaps();
this->optimizeToMeanDisplacement(mean_displacement); //espace admissible
this->computePressures(mean_displacement); //p=K-1*u
this->computeGaps();
G = this->computeG();
std::cout << "G vaut "<< G<< std::endl;
}
f = computeStoppingCriterion();
if (nb_iterations % dump_freq == 0) {
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << std::fixed << std::endl;
}
convergence_iterations.push_back(f);
// if(nb_iterations == 2) return 0;
}
this->computePressures(mean_displacement);
return f;
}
/* -------------------------------------------------------------------------- */
Real BemGigi::computeStoppingCriterion() {
Real crit=0.;
Real disp_norm = 0.;
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for reduction(+:crit, disp_norm)
for (UInt i = 0; i < size; ++i) {
crit += ( this->old_displacements(i) - true_displacements(i))*(
this->old_displacements(i) - true_displacements(i));
disp_norm += (true_displacements(i)*true_displacements(i));
}
return crit / disp_norm;
}
/* -------------------------------------------------------------------------- */
Real BemGigi::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->surface_tractions(i);
res += val*val;}
}
STOPTIMER("computeG");
return res;
}
/* -------------------------------------------------------------------------- */
void BemGigi::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 BemGigi::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 BemGigi::updateSearchDirection(Real factor) {
STARTTIMER("updateSearchDirection");
UInt n = surface.size();
UInt size = n*n;
Surface<Real> & gradF = this->functional->getGradF();
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
// this->search_direction(i) *=factor;
this->search_direction(i) = gradF(i);
}
STOPTIMER("updateSearchDirection");
}
/* -------------------------------------------------------------------------- */
Real BemGigi::computeOptimalStep() {
STARTTIMER("computeOptimalStep");
search_direction.FFTTransform(spectral_search_direction, nthreads);
UInt n = surface.size();
UInt size = n*n;
spectral_search_direction(0) = 0;
#pragma omp parallel for
for (UInt i = 1; i < size; ++i) {
spectral_search_direction(i) /= this->surface_spectral_influence_disp(i);
}
// /!\ does not contain the spectral search direction anymore
spectral_search_direction.FFTITransform(nthreads);
// Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
// spectral_search_direction -= average;
Surface<Real> & gradF = this->functional->getGradF();
Real numerator = 0., denominator = 0.;
#pragma omp parallel for reduction(+: numerator, denominator)
for (UInt i = 0; i < size; ++i) {
numerator += gradF(i) * search_direction(i);
denominator += spectral_search_direction(i).real() * search_direction(i);
}
Real alpha = numerator / denominator;
STOPTIMER("computeOptimalStep");
return alpha;
}
/* -------------------------------------------------------------------------- */
Real BemGigi::computeTau() {
STARTTIMER("computeOptimalStep");
search_direction.FFTTransform(spectral_search_direction, nthreads);
UInt n = surface.size();
UInt size = n*n;
spectral_search_direction(0) = 0;
#pragma omp parallel for
for (UInt i = 1; i < size; ++i) {
spectral_search_direction(i) /= this->surface_spectral_influence_disp(i);
}
// /!\ does not contain the spectral search direction anymore
spectral_search_direction.FFTITransform(nthreads);
// Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
// spectral_search_direction -= average;
Surface<Real> & gradF = this->functional->getGradF();
Real numerator = 0., denominator = 0.;
#pragma omp parallel for reduction(+: numerator, denominator)
for (UInt i = 0; i < size; ++i) {
if (this->gap(i) > 0) {
numerator += gradF(i) * search_direction(i);
denominator += spectral_search_direction(i).real() * search_direction(i);}
}
Real tau = numerator / denominator;
STOPTIMER("computeOptimalStep");
return tau;
}
/* -------------------------------------------------------------------------- */
void BemGigi::updateDisplacements(Real alpha) {
STARTTIMER("updateDisplacements");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->true_displacements(i) -= alpha*this->search_direction(i);
}
STOPTIMER("updateDisplacements");
}
/* -------------------------------------------------------------------------- */
void BemGigi::updateDisplacementsB(Real alpha) {
STARTTIMER("updateDisplacements");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if (this->gap(i)>0){
this->true_displacements(i) -= alpha*this->search_direction(i);
}
}
STOPTIMER("updateDisplacements");
}
/* -------------------------------------------------------------------------- */
void BemGigi::emptyOverlap() {
STARTTIMER("emptyoverlap");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if (gap(i) < 0)
this->true_displacements(i) = this->surface(i);
}
STOPTIMER("emptyoverlap");
}
/* -------------------------------------------------------------------------- */
void BemGigi::enforceMeanDisplacement(Real mean_displacement) {
STARTTIMER("enforceMeanDisplacement");
UInt n = surface.size();
UInt size = n*n;
Real average = SurfaceStatistics::computeAverage(this->true_displacements, 0);
Real factor = mean_displacement / average;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->true_displacements(i) *= factor;
}
STOPTIMER("enforceMeanDisplacement");
}
/* -------------------------------------------------------------------------- */
void BemGigi::computePressures(Real mean_displacement) {
UInt n = surface.size();
UInt size = n*n;
Real moy_surface=SurfaceStatistics::computeAverage(this->surface, 0);
this->computeTractionsFromDisplacements();
Real true_pressure = 0.;
for (UInt i = 0 ; i < n * n ; i++) {
true_pressure+=this->surface_tractions(i)*(this->true_displacements(i)-mean_displacement-this->surface(i));
}
true_pressure/=(moy_surface*size-mean_displacement*size);
std::cout << "true_pressure "<< true_pressure << std::endl;
//this->surface_tractions += true_pressure;
}
/* -------------------------------------------------------------------------- */
void BemGigi::computeTruePressures(Real mean_displacement) {
this->computeGaps();
this->functional->computeGradF();
UInt n = gap.size();
//Orthogonalite
for (UInt i = 0 ; i < n * n ; i++) {
if (this->gap(i) > 0){
this->surface_tractions(i)=0;
}
}
}
Event Timeline
Log In to Comment