diff --git a/src/bem_polonski.cpp b/src/bem_polonski.cpp index 3f3e0ab..c104009 100644 --- a/src/bem_polonski.cpp +++ b/src/bem_polonski.cpp @@ -1,300 +1,298 @@ /** * * @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_polonski.hh" #include "fft_plan_manager.hh" #include "fftransform.hh" #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define TIMER #include "surface_timer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ BemPolonski::BemPolonski(Surface & p) : BemFFTBase(p), surface_t(p.size(),p.getL()), surface_r(p.size(),p.getL()), - surface_spectral_r(p.size(),p.getL()), pold(p.size(),p.getL()){ e_star = 1.; max_iterations = 10000; } /* -------------------------------------------------------------------------- */ BemPolonski::~BemPolonski() {} /* -------------------------------------------------------------------------- */ Real BemPolonski::computeEquilibrium(Real epsilon, Real pressure) { this->computeSpectralInfluenceOverDisplacement(); this->surface_t = 0.; this->surface_r = 0.; - this->surface_spectral_r = 0.; this->pold = 0.; this->surface_displacements = 0.; Real delta = 0.; Real Gold = 1.; Real f = 1e300; Real fPrevious = 1e300; Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions); if (current_pressure <= 0.) surface_tractions = pressure; this->enforcePressureBalance(pressure); convergence_iterations.clear(); nb_iterations = 0; while(f > epsilon && nb_iterations < max_iterations) { fPrevious = f; this->computeDisplacementsFromTractions(); this->computeGaps(); Real gbar = this->computeMeanGapsInContact(); this->gap -= gbar; Real G = this->computeG(); //std::cout << "G vaut "<< G<< std::endl; this->updateT(G,Gold,delta); Real tau = this->computeTau(); pold = this->surface_tractions; Gold = G; delta = this->updateTractions(tau); //Projection on admissible space this->enforcePressureBalance(pressure); f = this->computeF(); if (nb_iterations % dump_freq == 0){ Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions); std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f-fPrevious << " " << A << std::fixed << std::endl; } convergence_iterations.push_back(f); ++nb_iterations; } return f; } /* -------------------------------------------------------------------------- */ void BemPolonski::computeGaps() { UInt size = surface.size(); #pragma omp parallel for for (UInt i=0; i < size*size; i++) { gap(i) = surface_displacements(i) - surface(i); } } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeMeanGapsInContact(){ STARTTIMER("computeMeanGapsInContact"); UInt n = surface.size(); UInt size = n*n; Real res = 0.; UInt nb_contact = 0; #pragma omp parallel for reduction(+: nb_contact,res) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; ++nb_contact; res += this->gap(i); } res /= nb_contact; STOPTIMER("computeMeanGapsInContact"); return res; } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeG(){ STARTTIMER("computeG"); UInt n = surface.size(); UInt size = n*n; Real res = 0.; #pragma omp parallel for reduction(+: res) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; Real val = this->gap(i); res += val*val; } STOPTIMER("computeG"); return res; } /* -------------------------------------------------------------------------- */ void BemPolonski::updateT(Real G, Real Gold, Real delta){ STARTTIMER("updateT"); UInt n = surface.size(); UInt size = n*n; Real factor = delta*G/Gold; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) this->surface_t(i) = 0.; else{ this->surface_t(i) *= factor; this->surface_t(i) += this->gap(i); } } STOPTIMER("updateT"); } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeTau(){ STARTTIMER("computeTau"); this->applyInfluenceFunctions(surface_t, surface_r); const UInt size = surface_t.size() * surface_t.size(); Real rbar = 0; UInt nb_contact = 0; #pragma omp parallel for reduction(+: nb_contact,rbar) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(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 (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; tau_sum1 += this->gap(i)*surface_t(i); tau_sum2 += surface_r(i)*surface_t(i); } Real tau = tau_sum1/tau_sum2; STOPTIMER("computeTau"); return tau; } /* -------------------------------------------------------------------------- */ Real BemPolonski::updateTractions(Real tau){ STARTTIMER("updateTractions"); UInt n = surface.size(); UInt size = n*n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; this->surface_tractions(i) -= tau*this->surface_t(i); if (this->surface_tractions(i) < 0.) this->surface_tractions(i) = 0; } //compute number of interpenetration without contact UInt nb_iol = 0; #pragma omp parallel for reduction(+: nb_iol) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0. && this->gap(i) < 0.) { this->surface_tractions(i) -= tau*this->gap(i); ++nb_iol; } } Real delta = 0; if (nb_iol > 0) delta = 0.; else delta = 1.; STOPTIMER("updateTractions"); return delta; } /* -------------------------------------------------------------------------- */ void BemPolonski::enforcePressureBalance(Real applied_pressure){ STARTTIMER("enforcePressureBalance"); UInt n = surface.size(); UInt size = n*n; Real pressure = 0.; #pragma omp parallel for reduction(+: pressure) for (UInt i = 0; i < size; ++i) { pressure += this->surface_tractions(i); } pressure *= 1./size; this->surface_tractions *= applied_pressure/pressure; STOPTIMER("enforcePressureBalance"); } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeF() { UInt n = surface.size(); UInt size = n*n; Real res = 0; Real d_max = SurfaceStatistics::computeMaximum(surface_displacements); Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions)); Real s_max = SurfaceStatistics::computeMaximum(surface); #pragma omp parallel for reduction(+:res) for (UInt i = 0; i < size; ++i) { res += std::abs(surface_tractions(i) * (surface_displacements(i)-d_max-surface(i)+s_max)); // res += // this->surface_tractions[i].real() // *(surface_displacements[i].real() - surface[i].real()); } return res/t_sum; } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/bem_polonski.hh b/src/bem_polonski.hh index c276272..2b752bc 100644 --- a/src/bem_polonski.hh +++ b/src/bem_polonski.hh @@ -1,99 +1,97 @@ /** * * @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_POLONSKI_H #define BEM_POLONSKI_H /* -------------------------------------------------------------------------- */ #include "bem_fft_base.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ class BemPolonski : public BemFFTBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: BemPolonski(Surface & p); virtual ~BemPolonski(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: //! compute the equilibrium situation Real computeEquilibrium(Real epsilon,Real pressure); //! compute the mean of the gaps in the contact region Real computeMeanGapsInContact(); //! compute the non zero gap penalty functional Real computeG(); //! update the search direction void updateT(Real G, Real Gold, Real delta); //! compute the conjugate step size Real computeTau(); //! update the tractions Real updateTractions(Real tau); //! enforce the applied and contact pressures balance void enforcePressureBalance(Real applied_pressure); //! compute the gaps void computeGaps(); //! compute F Real computeF(); /* ------------------------------------------------------------------------ */ /* 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;}; Surface & getSurfaceT(){return surface_t;}; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! exploration direction for the polonski algo Surface surface_t; //! projected exploration direction for the polonski algo Surface surface_r; - //! projected exploration direction for the polonski algo (fourier version) - SurfaceComplex surface_spectral_r; //! previous pressure Surface pold; }; __END_TAMAAS__ #endif /* BEM_POLONSKI_H */