Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/bem/bem_gigipol.cpp
a049dab2b91bfeatures/netcdf-mpi

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 "bem_gigipol.hh"
#include "surface.hh"
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <vector>
/* -------------------------------------------------------------------------- */
__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() {
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;
}
}
return res;
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
Real BemGigipol::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;
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) {
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);
}
}
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeTau() {
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;
return tau;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeTau(Real mean_pressure) {
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;
return tau;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::updateDisplacements(Real tau, Real mean_displacement) {
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;
}
/* -------------------------------------------------------------------------- */
void BemGigipol::enforceMeanDisplacement(Real mean_displacement) {
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);
}
}
/* -------------------------------------------------------------------------- */
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