Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F64408941
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 26, 16:44
Size
8 KB
Mime Type
text/x-c
Expires
Tue, May 28, 16:44 (2 d)
Engine
blob
Format
Raw Data
Handle
17898560
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"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
Real BemGigi::computeEquilibrium(Real epsilon,
Real mean_displacement) {
this->computeSpectralInfluenceOverDisplacement();
Real Rold = 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();
convergence_iterations.clear();
nb_iterations = 0;
max_iterations =50;
Real R=0;
while (f> epsilon && nb_iterations++ < max_iterations) {
this->computeGaps();
this->functional->computeGradFU();
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);
this->computeGaps();
this->optimizeToMeanDisplacement(mean_displacement); //espace admissible
this->computeGaps();
this->computePressures(mean_displacement);
f=computeStoppingCriterion();
}
this->computePressures(mean_displacement);
return f;
}
/* -------------------------------------------------------------------------- */
Real BemGigi::computeStoppingCriterion() {
Real rho = this->functional->getParameter("rho");
Real surface_energy = this->functional->getParameter("surface_energy");
Real crit=0.;
//Real disp_norm = 0.;
UInt n = surface.size();
UInt size = n*n;
//#pragma omp parallel for reduction(+:crit, disp_norm)
#pragma omp parallel for reduction(+:crit)
for (UInt i = 0; i < size; ++i) {
crit += std::abs(this->gap(i)*(this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho)));
}
return crit;
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
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;
const Surface<Real> & gradF = this->functional->getGradF();
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->search_direction(i)= gradF(i);
}
STOPTIMER("updateSearchDirection");
}
/* -------------------------------------------------------------------------- */
Real BemGigi::computeOptimalStep() {
STARTTIMER("computeOptimalStep");
this->applyInverseInfluenceFunctions(search_direction, projected_direction);
UInt n = surface.size();
UInt size = n*n;
const 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 += projected_direction(i) * search_direction(i);
}
Real alpha = numerator / denominator;
STOPTIMER("computeOptimalStep");
return alpha;
}
/* -------------------------------------------------------------------------- */
Real BemGigi::computeTau() {
return computeOptimalStep();
}
/* -------------------------------------------------------------------------- */
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::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 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 BemGigi::computePressures(Real mean_displacement) {
this->computeTractionsFromDisplacements();
UInt n = surface.size();
UInt size = n*n;
Real rho = this->functional->getParameter("rho");
Real surface_energy = this->functional->getParameter("surface_energy");
Real mini=3000;
for (UInt i = 0; i < size; ++i) {
if (gap(i)<rho){
if (this->surface_tractions(i)+surface_energy/rho <mini )
mini= this->surface_tractions(i)+surface_energy/rho ;
}
else
{ if (this->surface_tractions(i)<mini )
mini= this->surface_tractions(i) ;}
}
this->surface_tractions-=mini;
}
__END_TAMAAS__
Event Timeline
Log In to Comment