diff --git a/src/bem_kato.cpp b/src/bem_kato.cpp index 5133229..76c2cd8 100644 --- a/src/bem_kato.cpp +++ b/src/bem_kato.cpp @@ -1,253 +1,273 @@ /** * * @author Guillaume Anciaux * * @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 . * */ /* -------------------------------------------------------------------------- */ #include #include "surface.hh" #include "bem_kato.hh" #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ Real BemKato::linescan(Real scale,Real pressure){ updatePressure(scale); shiftPressure(pressure); truncatePressure(); this->functional->computeF(); Real res = this->functional->getF(); return res; } /* -------------------------------------------------------------------------- */ Real BemKato::linesearch(Real & hmax, Real fold, Real pressure, int search_flag){ if (search_flag){ Real h = hmax; // Real fold = bem.computeF(); //if (fold == 0) fold = 1e300; Real f = linescan(h,pressure); if (f < fold) return f; while (f > fold){ h *= 0.5; if (h < 1e-3) throw 1; f = linescan(h,pressure); } f = linescan(h,pressure); // if (hmax / h > 10) hmax /=2; return f; } return linescan(hmax,pressure); } /* -------------------------------------------------------------------------- */ Real BemKato::computeEquilibrium(Real epsilon, Real pressure){ // UInt n = surface.size(); // UInt size = n*n; + this->computeSpectralInfluenceOverDisplacement(); this->computeDisplacementsFromTractions(); this->functional->computeGradF(); this->backupTractions(); Real f = 1.; Real fPrevious = 1e300; this->nb_iterations = 0; this->convergence_iterations.clear(); while(f > epsilon && this->nb_iterations < this->max_iterations) { fPrevious = f; this->computeDisplacementsFromTractions(); this->functional->computeGradF(); this->backupTractions(); try { f = linescan(1.,pressure); } catch (int e){ std::cout << " line search failed " << std::endl; f = linescan(1,pressure); nb_iterations = 0; break; } if (nb_iterations % dump_freq == 0) { // std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f-fold << " " << ((f-fold)/forigin) << std::endl; std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f-fPrevious << std::endl; } convergence_iterations.push_back(f); ++nb_iterations; } this->computeTrueDisplacements(); this->computeGaps(); return f; } /* -------------------------------------------------------------------------- */ void BemKato::updatePressure(Real scale){ STARTTIMER("updatePressure"); UInt n = surface.size(); UInt size = n*n; const Surface & gradF = this->functional->getGradF(); #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->surface_tractions(i) = this->surface_tractions_backup(i) - gradF(i) * scale; } STOPTIMER("updatePressure"); } /* -------------------------------------------------------------------------- */ void BemKato::backupTractions(){ STARTTIMER("switchPressure"); UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->surface_tractions_backup(i) = this->surface_tractions(i); } STOPTIMER("switchPressure"); } /* -------------------------------------------------------------------------- */ Real BemKato::positivePressure(Real step){ STARTTIMER("positivePressure"); UInt n = surface.size(); UInt size = n*n; Real p_tot = 0.0; #pragma omp parallel for reduction(+:p_tot) for (UInt i = 0; i < size; ++i) { Real sh_press = this->surface_tractions(i) + step; if (sh_press > max_pressure) p_tot += max_pressure; else p_tot += sh_press*(sh_press > 0); } STOPTIMER("positivePressure"); return p_tot/n/n; } /* -------------------------------------------------------------------------- */ void BemKato::shiftPressure(Real required_pressure){ Real step_min = -10; Real step_max = 10; STARTTIMER("shiftPressureInitialGuess"); Real p_max = positivePressure(step_max); Real p_min = positivePressure(step_min); for(UInt i = 0; p_max < required_pressure && i < 8; ++i, step_max*=10) { p_max = positivePressure(step_max); } for(UInt i = 0; p_min > required_pressure && i < 8; ++i, step_min*=10) { p_min = positivePressure(step_min); } Real p = positivePressure(0.0); Real epsilon = 1e-12; STOPTIMER("shiftPressureInitialGuess"); STARTTIMER("shiftPressureDichotomy"); while (fabs(step_min - step_max) > epsilon){ Real step = (step_min + step_max)/2; p = positivePressure(step); if (p > required_pressure) step_max = step; else if (p < required_pressure) step_min = step; else { step_max = step; step_min = step; } } STOPTIMER("shiftPressureDichotomy"); STARTTIMER("shiftPressureTrueShift"); UInt n = surface.size(); UInt size = n*n; //shift the pressure so that satisfies the constraint #pragma omp parallel for for (UInt i = 0; i < size; ++i) { this->surface_tractions(i) += (step_max+step_min)/2; } STOPTIMER("shiftPressureTrueShift"); } /* -------------------------------------------------------------------------- */ void BemKato::truncatePressure(){ STARTTIMER("truncatePressure"); UInt n = surface.size(); UInt size = n*n; //shift the pressure so that satisfies the constraint #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) > max_pressure) this->surface_tractions(i) = max_pressure; else this->surface_tractions(i) *= (this->surface_tractions(i) > 0); } STOPTIMER("truncatePressure"); } /* -------------------------------------------------------------------------- */ +void BemKato::computeTrueDisplacements() { + UInt n = surface.size(); + UInt size = n*n; + + Real shift = 1e300; +#pragma omp parallel for reduce(min: shift) + for (UInt i = 0 ; i < size ; i++) { + if (surface_displacements(i) - shift - surface(i) < 0. && + surface_tractions(i) < max_pressure){ + shift = surface_displacements(i) - surface(i); + } + } + + this->true_displacements = surface_displacements; + this->true_displacements -= shift; +} + +/* -------------------------------------------------------------------------- */ + __END_TAMAAS__ diff --git a/src/bem_kato.hh b/src/bem_kato.hh index 1a48985..63cd877 100644 --- a/src/bem_kato.hh +++ b/src/bem_kato.hh @@ -1,104 +1,106 @@ /** * * @author Guillaume Anciaux * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef BEM_KATO_H #define BEM_KATO_H /* -------------------------------------------------------------------------- */ #include "bem_fft_base.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class BemKato : public BemFFTBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: BemKato(Surface & p) : BemFFTBase(p), surface_tractions_backup(p.size(),p.getL()){ e_star = 1.; max_pressure = 1e300; }; virtual ~BemKato(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: //! compute the equilibrium situation Real computeEquilibrium(Real epsilon,Real pressure); public: //! perform a line scan in a given direction Real linescan(Real scale,Real pressure); //! perform a line search Real linesearch(Real & hmax, Real fold, Real pressure, int search_flag); //! make a step in the gradient direction(copy in pressure2) void updatePressure(Real scale); //! backup tractions before doing the linescan void backupTractions(); //! find the constant pressure shift making balance with load void shiftPressure(Real required_pressure); //! truncate the pressure to positive values Real positivePressure(Real step); //! compute the reaction (sum positive pressures) void truncatePressure(); + //! compute true displacements + void computeTrueDisplacements(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: //! get the number of iterations UInt getNbIterations() const{return this->nb_iterations;}; //! get the convergence steps of iterations const std::vector & getConvergenceIterations() const{return this->convergence_iterations;}; //! set the maximal pressure admissible void setMaxPressure(Real max_pressure){this->max_pressure = max_pressure;}; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! surface of tractions used during steepest descent Surface surface_tractions_backup; //! maximal pressure admissible (for perfectly plastic purpose) Real max_pressure; }; __END_TAMAAS__ #endif /* BEM_KATO_H */