f[x_] := Exp[2/(Abs[x - 1] + 1)]*Exp[-(x/2)^6]; (* artificially produced cusped function; in actual case this will produced numerically *)

xi = -4; (* initial point on the grid *)
xf = +4; (* final point on the grid *)
nPts = 1000; (* nr of points on the grid *)
d = (xf - xi)/(nPts - 1); (* distance between two points on the (uniform) grid *)
data = Table[N[{xi + i*d, f[a + i*d]}], {i, 0, nPts - 1}]; (* grid points *)
splFunct = Interpolation[data, Method -> {"Spline"}]; (* interpolating spline from the numerically produced data points *)
g[x_] := Piecewise[{{splFunct[x], xi <= x <= xf}}, 0]; (* a function which is equivalent to the spline function on the interval [xi , xf], and vanishes outside *)

intg[x_, c1_, c2_] := g[x]*g[x - c1]*g[x - 2 c2]; (* integrand made from a product of functions with different offsets in the argument, c1 and c2 *)

(* integration for each point on the grid; the result is stored in 'intgResult': *)
Timing[
    intgResult = Table[{data[[i, 1]],
     NIntegrate[intg[x, data[[i, 1]], data[[i, 1]]], {x, xi, xf},
      Method -> {Automatic, "SymbolicProcessing" -> 0},
      AccuracyGoal -> 10, MaxPoints -> 2000]},
      {i, 1, Length[data]}];
    ]