diff --git a/src/model/influence.hh b/src/model/influence.hh
index 40e8f7d..dad2149 100644
--- a/src/model/influence.hh
+++ b/src/model/influence.hh
@@ -1,299 +1,303 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 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 "loop.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 namespace influence {
 
 /// Mother class for influence functions
 template <UInt dim>
 class Influence {
 public:
   Influence(VectorProxy<const Real, dim> domain) : domain(domain) {}
 
 protected:
   VectorProxy<const Real, dim> domain;
 };
 
 /// Class representing the basic influence functions from Stanley & Kato (1997)
 template <UInt dim>
 class Basic : public Influence<dim> {
 public:
   Basic(Real E_hertz, VectorProxy<const Real, dim> domain)
       : Influence<dim>(domain), E_hertz(E_hertz) {}
 
   /// Influence function: computes the influence coefficient F
   CUDA_LAMBDA void operator()(VectorProxy<Real, dim>&& q_vec,
                               Complex& F) const {
     q_vec *= 2 * M_PI;
     q_vec /= this->domain;
     F = 2. / (q_vec.l2norm() * E_hertz);
   }
 
 private:
   Real E_hertz;
 };
 
 template <UInt dim>
 class Surface : public Influence<dim> {};
 
 /// Class for the influence with all components of applied tractions (dim = 2)
 template <>
 class Surface<2> : public Influence<2> {
   constexpr static UInt dim = 2;
 
 public:
   Surface(Real E, Real nu, VectorProxy<const Real, dim> domain)
       : Influence<dim>(domain), E(E), nu(nu) {}
 
   CUDA_LAMBDA void
   operator()(VectorProxy<Real, dim>&& q_,
              MatrixProxy<Complex, dim + 1, dim + 1>&& F) const {
     q_ *= 2 * M_PI;
     q_ /= this->domain;
     const Real q = q_.l2norm();
     q_ /= q;
 
     // Diagonal
     F(0, 0) = 2 * (1 + nu) * (1 - nu * q_(0) * q_(0));
     F(1, 1) = 2 * (1 + nu) * (1 - nu * q_(1) * q_(1));
     F(2, 2) = 2 * (1 - nu * nu);
 
     F(0, 1) = F(1, 0) = -q_(0) * q_(1) * 2 * nu * (1 + nu);
 
     F(0, 2) = I * q_(0) * (1 + nu) * (1 - 2 * nu);
     F(1, 2) = I * q_(1) * (1 + nu) * (1 - 2 * nu);
 
     F(2, 0) = -F(0, 2);
     F(2, 1) = -F(1, 2);
 
     F *= 1. / (E * q);
   }
 
 private:
   Real E, nu;
   const Complex I = Complex(0, 1);
 };
 
 /// Class for the influence with all components of applied tractions (dim = 2)
 template <>
 class Surface<1> : public Influence<1> {
   constexpr static UInt dim = 1;
 
 public:
   Surface(Real E, Real nu, VectorProxy<const Real, dim> domain)
       : Influence<dim>(domain), E(E), nu(nu) {}
 
   CUDA_LAMBDA void
   operator()(VectorProxy<Real, dim>&& q_,
              MatrixProxy<Complex, dim + 1, dim + 1>&& F) const {
     q_ *= 2 * M_PI;
     q_ /= this->domain;
     const Real q = q_.l2norm();
 
     // First line of influence matrix
     F(0, 0) = 2. / q * (1 - nu * nu);
     F(0, 1) = I / q_(0) * (1 + nu) * (1 - 2 * nu);
 
     // Second line of influence matrix
     F(1, 0) = -F(0, 1);
     F(1, 1) = F(0, 0);
 
     F *= 1. / E;
   }
 
 private:
   Real E, nu;
   const Complex I = Complex(0, 1);
 };
 
 /// Class for the Kelvin tensor
 template <UInt dim, UInt derivative_order>
 class Kelvin;
 
 template <UInt derivative_order>
 class Kelvin<3, derivative_order> : Influence<3> {
   constexpr static UInt dim = 3;
 
 public:
   Kelvin(Real mu, Real nu, VectorProxy<const Real, dim> domain)
       : Influence<dim>(domain), mu(mu), nu(nu) {}
 
 private:
   /// Compute leading linear term A
   template <bool upper>
   inline void computeA(VectorProxy<const Real, dim - 1>& q_,
                        Matrix<Complex, dim, dim>& A) const;
   /// Compute leading constant term B, which does not depend on y_3
   inline void computeB(VectorProxy<const Real, dim - 1>& q_,
                        Matrix<Complex, dim, dim>& B) const;
 
 public:
   /// Compute leading terms
   template <bool upper>
   inline void computeTerms(VectorProxy<const Real, dim - 1>& q_,
                            Matrix<Complex, dim, dim>& A,
                            Matrix<Complex, dim, dim>& B) const {
     computeA<upper>(q_, A);
     computeB(q_, B);
   }
 
 private:
   Real mu, nu;
   const Complex I = Complex(0, 1);
 };
 
 /* -------------------------------------------------------------------------- */
 /* Base Kelvin tensor */
 /* -------------------------------------------------------------------------- */
 
 template <>
 inline void Kelvin<3, 0>::computeB(VectorProxy<const Real, dim - 1>& q_,
                                    Matrix<Complex, dim, dim>& B) const {
   const Real q = q_.l2norm();
   const Real bq1 = q_(0) / q;
   const Real bq2 = q_(1) / q;
 
   for (UInt i = 0; i < dim; ++i)
     B(i, i) = 4 * (1 - nu);
   B(2, 2) -= 1;
 
   B(0, 0) = -bq1 * bq1;
   B(1, 1) = -bq2 * bq2;
   B(0, 1) = B(1, 0) = -bq1 * bq2;
 
   B(0, 2) = B(1, 2) = B(2, 0) = B(2, 1) = 0;
   B *= (1 / (8 * mu * q * (1 - nu)));
 }
 
 template <>
 template <bool upper>
 inline void Kelvin<3, 0>::computeA(VectorProxy<const Real, dim - 1>& q_,
                                    Matrix<Complex, dim, dim>& A) const {
   const Real q = q_.l2norm();
   const Real bq1 = q_(0) / q;
   const Real bq2 = q_(1) / q;
 
   constexpr Real sign = (upper) ? 1. : -1.;
 
   A(0, 0) = sign * bq1 * bq1;
   A(1, 1) = sign * bq2 * bq2;
   A(2, 2) = sign;
 
   A(1, 0) = A(0, 1) = sign * bq1 * bq2;
   A(2, 0) = A(0, 2) = -I * bq1;
   A(2, 1) = A(1, 2) = -I * bq2;
 
   A *= (1 / (8 * mu * (1 - nu)));
 }
 
 /* -------------------------------------------------------------------------- */
 
 class KelvinIntegrator {
   static inline Complex constantPrimitive(const Real& x, const Real& alpha,
                                           const Complex& beta,
                                           const Complex& gamma) {
     return std::exp(alpha * x) / alpha * (beta * x + gamma - beta / alpha);
   }
 
   static inline Complex linearPrimitive(const Real& x, const Real& alpha,
                                         const Complex& beta,
                                         const Complex& gamma) {
     auto coeff = gamma - 2. * beta / alpha;
     return std::exp(alpha * x) / alpha *
            (beta * x * x + coeff * x - coeff / alpha);
   }
 
   static inline void computePolyCoefficients(Complex& beta, Complex& gamma,
                                              const Real& dl, const Real& dij,
                                              const Complex& A,
                                              const Complex& B) {
     beta = A * dl;
     gamma = A * dij + B;
   }
 
 public:
   template <bool upper, UInt dim>
   static inline void partialIntegration(Matrix<Complex, dim, dim>& F,
                                         const Real& alpha, const Real& dij,
                                         const Real& dl,
                                         const Matrix<Complex, dim, dim>& A,
                                         const Matrix<Complex, dim, dim>& B) {
     Real xsup = (upper) ? 1 : 0;
     Real xinf = (upper) ? 0 : -1;
     Real multiplier = (upper) ? -1 : 1;
     Real e_t = std::exp(alpha * dij / dl);
     Complex beta = 0, gamma = 0;
     for (UInt i = 0; i < dim; ++i) {
       for (UInt j = 0; j < dim; ++j) {
         computePolyCoefficients(beta, gamma, dl, dij, A(i, j), B(i, j));
         F(i, j) += e_t * (constantPrimitive(xsup, alpha, beta, gamma) -
                           constantPrimitive(xinf, alpha, beta, gamma));
         F(i, j) +=
             e_t * multiplier * (linearPrimitive(xsup, alpha, beta, gamma) -
                                 linearPrimitive(xinf, alpha, beta, gamma));
       }
     }
   }
 
   template <UInt dim, UInt derivative, Int yj_xi>
   static inline void integrate(const Kelvin<dim, derivative>& kelvin, Real yj,
                                Real xi, Real dl,
                                VectorProxy<const Real, dim - 1>& q,
                                Matrix<Complex, dim, dim>& F) {
     F = 0;
     Matrix<Complex, dim, dim> A, B;
 
     const Real q_norm = q.l2norm();
     const Real dij = yj - xi;
+
+    if (-q_norm * std::abs(dij) < std::log(1e-2))
+      return;
+
     constexpr Real sign = (yj_xi > 0) ? 1. : -1;
 
     Real alpha = 0;
 
     // Integrate [-1, 0]
     kelvin.template computeTerms<(yj_xi > 0)>(q, A, B);
     alpha = -sign * q_norm * dl;
 
     partialIntegration<false>(F, alpha, dij, dl, A, B);
 
     // Integrate [0, 1]
     if (yj_xi == 0) {  // avoid computing the same A & B twice
       kelvin.template computeTerms<true>(q, A, B);
       alpha *= -1;
     }
 
     partialIntegration<true>(F, alpha, dij, dl, A, B);
   }
 };
 
 }  // namespace influence
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__