\documentclass{article}

\usepackage{amsmath}
\newcommand{\te}{\mathrm{e}}
\newcommand{\ti}{\mathrm{i}}
\newcommand{\td}{\mathrm{d}}

\begin{document}

\begin{equation}
g({\bf x},{\bf y}) = \frac{\te^{-\ti k r}}{4\pi r}, \qquad
r = |{\bf x} - {\bf y}|
\end{equation}

\begin{equation}
\frac{\partial^2 g({\bf x},{\bf y})}{\partial n_{\bf y}\partial n_{\bf x}} =
\nabla_{\bf y} \left( \nabla_{\bf x} g({\bf x},{\bf y}) {\bf n}_{\bf x} \right) {\bf n}_{\bf y}
\end{equation}

from Fischer and Matlab we know that:
%
\begin{equation}
\nabla_{\bf y} \left( \nabla_{\bf x} g({\bf x},{\bf y}) {\bf n}_{\bf x} \right) {\bf n}_{\bf y}
=
({\bf n}_x {\bf n}_y) (\nabla_{\bf x}\cdot \nabla_{\bf y} g({\bf x},{\bf y}))
-
{\bf n}_y \cdot \left[\nabla_{\bf y} \times \left({\bf n}_x\times \nabla_{\bf x}g({\bf x}, {\bf y})\right)\right]
\end{equation}
%
Since $\nabla_{\bf x} g = -\nabla_{\bf y} g$, the first term simplifies to $-\nabla_{\bf x}^2 g = k^2 g$:
%
\begin{equation}
\nabla_{\bf y} \left( \nabla_{\bf x} g({\bf x},{\bf y}) {\bf n}_{\bf x} \right) {\bf n}_{\bf y}
=
k^2 ({\bf n}_x {\bf n}_y) g({\bf x},{\bf y})
-
{\bf n}_y \cdot \left[\nabla_{\bf y} \times \left({\bf n}_x\times \nabla_{\bf x}g({\bf x}, {\bf y})\right)\right]
\end{equation}


\begin{multline}
\int_{S_x} p({\bf x}) \int_{S_y}
q({\bf y})
\frac{\partial^2 g({\bf x},{\bf y})}{\partial n_{\bf y}\partial n_{\bf x}}
\td S_y \td S_x
\\
= k^2 \int_{S_x} p({\bf x}) \int_{S_y}
q({\bf y})
({\bf n}_x {\bf n}_y) g({\bf x},{\bf y})
\td S_y \td S_x
\\
-\int_{S_x} p({\bf x})
\int_{S_y}
q({\bf y})
{\bf n}_y \cdot \left[\nabla_{\bf y} \times \left({\bf n}_x\times \nabla_{\bf x}g({\bf x}, {\bf y})\right)\right]
\td S_y \td S_x
\end{multline}
%
The first term on the rhs is weakly singular and can be evaluated easily.
The internal integral over $S_y$ in the second term can be transformed using the identity $\nabla \times (\alpha {\bf v}) = \nabla \alpha \times {\bf v} + \alpha \nabla \times {\bf v}$ and Stoke's theorem as
%
\begin{multline}
\int_{S_y} q({\bf y})
{\bf n}_y \cdot \left[\nabla_{\bf x} \times \left({\bf n}_x\times \nabla_{\bf x}g({\bf x}, {\bf y})\right)\right]
\td S_y \\
= \int_{S_y} q({\bf y})
\left[\nabla_{\bf y} \times \left({\bf n}_x\times \nabla_{\bf x}g({\bf x}, {\bf y})\right)\right]
\td {\bf S}_y \\
= \int_{S_y} \nabla_{\bf y} \times \left[
q({\bf y}) \left({\bf n}_x\times \nabla_{\bf x}g({\bf x}, {\bf y})\right)\right]
\td {\bf S}_y 
-
\int_{S_y} \nabla_{\bf y} q({\bf y}) \times 
\left({\bf n}_x\times \nabla_{\bf x}g({\bf x}, {\bf y})\right)
\td {\bf S}_y
\\
= \int_{\partial S_y} 
{\bf n}_x\times \nabla_{\bf x}\left(g({\bf x}, {\bf y}) q({\bf y})\right)
\td {\bf y}
-
\int_{S_y} \nabla_{\bf y} q({\bf y}) \times 
\left({\bf n}_x\times \nabla_{\bf x}g({\bf x}, {\bf y})\right)
{\bf n}_y \td S_y
\\
= \int_{\partial S_y} 
 \nabla_{\bf x}\left(g({\bf x}, {\bf y}) q({\bf y})\right)
\times \td {\bf y} \cdot {\bf n}_x
-
\int_{S_y} \left({\bf n}_y  \times \nabla_{\bf y} q({\bf y}) \right) \cdot
\left({\bf n}_x\times \nabla_{\bf x}g({\bf x}, {\bf y})\right)
\td S_y
\\
= \nabla_{\bf x} \times \int_{\partial S_y} 
g({\bf x}, {\bf y}) q({\bf y})
\td {\bf y} \cdot {\bf n}_x
-
\int_{S_y} \nabla_{\bf x}g({\bf x}, {\bf y}) \times \left({\bf n}_y  \times \nabla_{\bf y} q({\bf y}) \right) \cdot
{\bf n}_x
\td S_y
\\
= \nabla_{\bf x} \times
\left(\int_{\partial S_y} 
g({\bf x}, {\bf y}) q({\bf y})
\td {\bf y} 
-
\int_{S_y} g({\bf x}, {\bf y}) \left({\bf n}_y  \times \nabla_{\bf y} q({\bf y}) \right)
\td S_y
\right)
\cdot {\bf n}_x
\end{multline}
%
Introducing the result into the surface integral over $S_x$ and aplpying Stokes' theorem results in
%
\begin{multline}
\int_{S_x}
p({\bf x})
\nabla_{\bf x} \times
\left(\int_{\partial S_y} 
g({\bf x}, {\bf y}) q({\bf y})
\td {\bf y} 
-
\int_{S_y} g({\bf x}, {\bf y}) \left({\bf n}_y  \times \nabla_{\bf y} q({\bf y}) \right)
\td S_y
\right)
\cdot {\bf n}_x
\td S_x
\\
=
\int_{\partial S_x}
p({\bf x})
\left(\int_{\partial S_y} 
g({\bf x}, {\bf y}) q({\bf y})
\td {\bf y} 
-
\int_{S_y} g({\bf x}, {\bf y}) \left({\bf n}_y  \times \nabla_{\bf y} q({\bf y}) \right)
\td S_y
\right)
\td {\bf x}
\\
-
\int_{S_x}
\nabla_{\bf x} p({\bf x})
\times
\left(\int_{\partial S_y} 
g({\bf x}, {\bf y}) q({\bf y})
\td {\bf y} 
-
\int_{S_y} g({\bf x}, {\bf y}) \left({\bf n}_y  \times \nabla_{\bf y} q({\bf y}) \right)
\td S_y
\right)
\cdot {\bf n}_x
\td S_x
\\
=
\int_{\partial S_x}
p({\bf x})
\int_{\partial S_y} 
g({\bf x}, {\bf y}) q({\bf y})
\td {\bf y} 
\td {\bf x}
\\
-
\int_{\partial S_x}
p({\bf x})
\int_{S_y} g({\bf x}, {\bf y}) \left({\bf n}_y  \times \nabla_{\bf y} q({\bf y}) \right)
\td S_y
\td {\bf x}
\\
-
\int_{S_x}
\left({\bf n}_x \times \nabla_{\bf x} p({\bf x})\right)
\int_{\partial S_y} 
g({\bf x}, {\bf y}) q({\bf y})
\td {\bf y} 
\td S_x
\\
+
\int_{S_x}
\left({\bf n}_x \times \nabla_{\bf x} p({\bf x})\right)
\int_{S_y} g({\bf x}, {\bf y}) \left({\bf n}_y  \times \nabla_{\bf y} q({\bf y}) \right)
\td S_y
\td S_x
\end{multline}

\section{Computing $\nabla_{\bf x} N_i$}

\begin{equation}
\begin{bmatrix}
x'_\xi & y'_\xi & z'_\xi \\
x'_\eta & y'_\eta & z'_\eta \\
n_x & n_y & n_z
\end{bmatrix}
\left\{\begin{matrix}
N'_{i,x} \\ N'_{i,y} \\ N'_{i,z}
\end{matrix}\right\}
=
\left\{\begin{matrix}
N'_{i,\xi} \\ N'_{i,\eta} \\ 0
\end{matrix}\right\}
\end{equation}
%
or in compact form
%
\begin{equation}
\begin{bmatrix}
{{\bf r}'_\xi}^{\mathrm{T}} \\
{{\bf r}'_\eta}^{\mathrm{T}} \\
{\bf n}^{\mathrm{T}}
\end{bmatrix}
\nabla_{\bf x} N_i
=
\left\{\begin{matrix}
\nabla_{\xi} N_i \\ 0
\end{matrix}\right\}
\end{equation}
%
where ${\bf n} = {\bf r}'_\xi \times {\bf r}'_\eta$
%
Inverting the matrix yields
%
\begin{align}
\nabla_{\bf x} N_i
&= \frac{({\bf r}'_\eta\times {\bf n})N'_{i,\xi} + ({\bf n}\times{\bf r}'_\xi)  N'_{i,\eta}}{\| {\bf n} \|^2} \\
&= \frac{{\bf n} \times (N'_{i,\eta}{\bf r}'_\xi - N'_{i,\xi} {\bf r}'_\eta)}{\| {\bf n} \|^2} \\
&= \frac{{\bf n} \times \sum_j (N'_{i,\eta}  L'_{j,\xi} - N'_{i,\xi} L'_{j,\eta}) {\bf r}_j }{\| {\bf n} \|^2} \\
&= \frac{{\bf n} \times \sum_j Q_{ij} {\bf r}_j }{\| {\bf n} \|^2}
\end{align}

For the linear triangle element
%
\begin{equation}
{\bf Q} = \begin{bmatrix}
0 & 1 & -1\\
-1 & 0 & 1\\
1 & -1 & 0
\end{bmatrix}
\end{equation}


For the linear quad element
%
\begin{multline}
{\bf Q} = \frac{1}{8}\begin{bmatrix}
0 & \eta - 1 & \xi - \eta & 1 - \xi\\
1 - \eta & 0 &  - \xi - 1 & \eta + \xi\\
\eta - \xi & \xi + 1 & 0 &  - \eta - 1\\
\xi - 1 &  - \eta - \xi & \eta + 1 & 0
\end{bmatrix}
{\bf Q}_0 = \frac{1}{8}\begin{bmatrix}
0 & - 1 & 0 & 1 \\
1 & 0 &  - 1 & 0\\
0 & 1 & 0 &  - 1\\
- 1 &  0 &  1 & 0
\end{bmatrix} \\
{\bf Q}_\xi = \frac{1}{8}\begin{bmatrix}
0 & 0 & 1 & -1 \\
0 & 0 &  - 1 & 1\\
- 1 & 1 & 0 & 0\\
1 &  - 1 & 0 & 0
\end{bmatrix}
{\bf Q}_\eta = \frac{1}{8}\begin{bmatrix}
0 & 1  & - 1 & 0\\
- 1 & 0 &  0 & 1 \\
1  & 0 & 0 &  - 1\\
0 &  - 1 & 1 & 0
\end{bmatrix}
\end{multline}

\end{document}