  We study solutions to $Lu=f$ in $\Omega\subset\mathbb R^n$, being $L$ the
generator of any, possibly non-symmetric, stable L\'evy process.
  On the one hand, we study the regularity of solutions to $Lu=f$ in $\Omega$,
$u=0$ in $\Omega^c$, in $C^{1,\alpha}$ domains~$\Omega$. We show that solutions
$u$ satisfy $u/d^\gamma\in C^{\varepsilon_\circ}\big(\overline\Omega\big)$,
where $d$ is the distance to $\partial\Omega$, and $\gamma=\gamma(L,\nu)$ is an
explicit exponent that depends on the Fourier symbol of operator $L$ and on the
unit normal $\nu$ to the boundary $\partial\Omega$.
  On the other hand, we establish new integration by parts identities in half
spaces for such operators. These new identities extend previous ones for the
fractional Laplacian, but the non-symmetric setting presents some new
interesting features.
  Finally, we generalize the integration by parts identities in half spaces to
the case of bounded $C^{1,\alpha}$ domains. We do it via a new efficient
approximation argument, which exploits the H\"older regularity of $u/d^\gamma$.
This new approximation argument is interesting, we believe, even in the case of
the fractional Laplacian.
