RBF-OGr
Authors: Sage Shaw
Fri Feb 8 20:06:58 2019

This article outlines the RBF Orthogonal Gradients method for solving approximating surface operators as described in [1].

Overview

Given a set of points $\{\vec{x}_i\}_{i=1}^n$ from a surface $\SS$, and an operator $\mathcal{L}$, we would like to approximate the coressponding surface operator $\mathcal{L}_\SS$ of a function $u$ as a weighted sum of the function values $\{u(\vec{x}_i)\}_{i=1}^n$. In effect, we want the finite difference weights.

The idea is to construct an RBF interpolant $s(\vec{x})$ such that it interpolates the points, and has the property that $\mathcal{L}s(\vec{x}_i) = \mathcal{L}_\SS s(\vec{x}_i)$ for all $i=1, 2, ..., n$.

Three versions of this method are outlined here. The high order method, and low order method are from [1]. The symmetric method, is also by Dr. Piret, but has yet to be published.

Low Order Method

High Order Method

Symmetric Method

Overview

The idea is to set up an orthogonal gradients interpolant as before using only points on the surface (no isosurface reconstruction). From [2] we are given the lemma that $$ \Delta_\SS f = \Delta f - \frac{\partial^2 f}{\partial \vec{n}^2} - \kappa \frac{\partial f}{\partial \vec{n}} $$ where $\kappa$ is the mean curvature. For notational convenience let $\mathcal{L}_\vec{x} = \frac{\partial f}{\partial \vec{n}} = \vec{n}\cdot \nabla_\vec{x}$, and $\mathcal{H}_\vec{x} = \frac{\partial^2 f}{\partial \vec{n}^2} = \vec{n}^TH_\vec{x}(\cdot)\vec{n}$. Here $H_\vec{x}(\phi)$ denotes the hessian with respect to $\vec{x}$ of $\phi$ and $\vec{n}$ is the outward oriented unit normal vector to the surface at $\vec{x}$. Thus if we require that at a stencil center $\vec{y}$ we have $\mathcal{H}_\vec{x}s(\vec{x}) \vert_{\vec{x}=\vec{y}} = 0$ and $\mathcal{L}_\vec{x}s(\vec{x}) \vert_{\vec{x}=\vec{y}} = 0$ we will have the key property in the Orthogonal Gradients method: $\Delta s = \Delta_\SS s$, at least at $\vec{y}$.

Using a similar method to the moment conditions, we will Let $$ s(\vec{x}) = \sum c_i \phi(||\vec{x} - \vec{x}_i||) + d \mathcal{L}_\vec{y} \phi(||\vec{x} - \vec{y}||)+ e \mathcal{H}_\vec{y} \phi(||\vec{x} - \vec{y}||) $$ where $\vec{y} = \vec{x}_1$ is the stencil center, and require that $s$ interpolates the surface, $\mathcal{H}_\vec{x}s(\vec{x}) \vert_{\vec{x}=\vec{y}} = 0$, and $\mathcal{L}_\vec{x}s(\vec{x}) \vert_{\vec{x}=\vec{y}} = 0$.

Then we calculate a closed form for $\Delta s(\vec{x})$.

Lemmas

These lemmas use the same notation as specified in the Overview section. Additionally we we adopt the notation that for a given radial basis function $\phi$ (or $\phi_0$), define $\phi_1(r) = \frac{1}{r}\frac{d}{dr}\phi$ and for $i>1$ define $\phi_{i+1}(r) = \frac{1}{r}\frac{d}{dr}\phi_i$. The expression $\vec{n}\cdot(\vec{x}-\vec{y})$ occurs frequently so we will let $p = \vec{n}\cdot(\vec{x}-\vec{y})$.

Lemma 1

Denote $\vec{x} = \begin{bmatrix} x & y & z \end{bmatrix}^T, \vec{y} = \begin{bmatrix} a & b & c \end{bmatrix}^T$, and $r = ||\vec{x} - \vec{y}||$. Then \begin{align*} \frac{dr}{dx} &= \frac{d}{dx} \big((x-a)^2 + (y-b)^2 + (z-c)^2 \big)^{1/2} \\ &= 2(x-a)\frac{1}{2} \big((x-a)^2 + (y-b)^2 + (z-c)^2 \big)^{-1/2} \\ &= (x-a)\frac{1}{r} \\ \frac{d}{dx}\phi_n(r) &= \phi_n'(r)\frac{dr}{dx} = (x-a) \phi_{n+1}(r) \\ &\text{and similarly} \\ \nabla_\vec{x} \phi_n(r) &= (\vec{x}-\vec{y}) \phi_{n+1}(r) \\ \mathcal{L}_\vec{x}\phi_n(r) &= \vec{n}\cdot(\vec{x}-\vec{y}) \phi_{n+1}(r) \\ &= p\phi_{n+1}(r).\\ \end{align*} Also since $\frac{d}{da} \phi(r) = -(x-a)\phi_1(r) = - \frac{d}{dx}\phi(r)$ it follows that $\mathcal{L}_\vec{y} = -\mathcal{L}_\vec{x}$.

$\blacksquare$

Lemma 2

For $n \in \NN$, $$ \mathcal{L}_\vec{x}p^{n+1} = (n+1)p^n. $$

First denote $\vec{x} = \begin{bmatrix}x_1 & x_2 & x_3 \end{bmatrix}^T$ and $\vec{n} = \begin{bmatrix}n_1 & n_2 & n_3\end{bmatrix}$. Then \begin{align*} \frac{\partial}{\partial x_1}\big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n+1} &= (n+1) \big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n} n_1 \\ \frac{\partial}{\partial x_2}\big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n+1} &= (n+1) \big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n} n_2 \\ \frac{\partial}{\partial x_3}\big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n+1} &= (n+1) \big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n} n_3 \\ \nabla_\vec{x}\big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n+1} &= (n+1) \big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n} \vec{n} \\ \mathcal{L}_\vec{x} \big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n+1} &= \vec{n} \cdot \nabla_\vec{x}\big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n+1} \\ &= (n+1) \big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n} \vec{n} \cdot \vec{n} \\ &= (n+1) \big(\vec{n}\cdot(\vec{x}-\vec{y})\big)^{n} \end{align*}

$\blacksquare$

Lemma 3

For a radial basis function $\phi$, we have that $\mathcal{H}_\vec{x} \phi(r) = \mathcal{L}_\vec{x}\mathcal{L}_\vec{x} \phi(r) = \phi_1(r) + p^2 \phi_2(r)$. Also that $\mathcal{H}_\vec{x}\phi = \mathcal{H}_\vec{y}\phi$.

From Lemma 1 we have that $\mathcal{L}_\vec{x} \phi(r) = p \phi_1(r)$. Then \begin{align*} \mathcal{L}_\vec{x}\mathcal{L}_\vec{x} \phi(r) &= \mathcal{L}_\vec{x} \big( p \phi_1(r) \big) \\ &= \phi_1(r)\mathcal{L}_\vec{x} p + p \mathcal{L}_\vec{x} \phi_1(r) \\ &= \phi_1(r) + p^2 \phi_2(r). \\ \end{align*} Next consider \begin{align*} \frac{\partial}{\partial x} \phi(r) &= (x-a)\phi_1(r) \\ \frac{\partial^2}{\partial x^2} \phi(r) &= \phi_1(r) + (x-a)^2\phi_2(r) \\ \frac{\partial^2}{\partial xy} \phi(r) &= (x-a)(y-b)\phi_2(r). \end{align*} These will be useful in the calculation of the Hessian below \begin{align*} \mathcal{H}_\vec{x} \phi(r) &= \vec{n}^T \begin{bmatrix} \frac{\partial^2\phi}{\partial x^2} & \frac{\partial^2\phi}{\partial xy} & \frac{\partial^2\phi}{\partial xz} \\ \frac{\partial^2\phi}{\partial xy} & \frac{\partial^2\phi}{\partial y^2} & \frac{\partial^2\phi}{\partial yz} \\ \frac{\partial^2\phi}{\partial xz} & \frac{\partial^2\phi}{\partial yy} & \frac{\partial^2\phi}{\partial z^2} \end{bmatrix} \vec{n} \\ &= n_x^2\frac{\partial^2\phi}{\partial x^2} + n_y^2\frac{\partial^2\phi}{\partial y^2} + n_z^2\frac{\partial^2\phi}{\partial z^2} + 2 n_x n_y\frac{\partial^2\phi}{\partial xy} + 2 n_x n_z\frac{\partial^2\phi}{\partial xz} + 2 n_y n_z\frac{\partial^2\phi}{\partial yz} \\ &= (n_x^2 + n_y^2 + n_z^2)\phi_1(r) + \big( n_x(x-a) + n_y(y-b) + n_z(z-c) \big)^2\phi_2(r) \\ &= \phi_1(r) + p^2\phi_2(r). \end{align*} And lastly we have that $\mathcal{H}_\vec{x}\phi = \mathcal{L}_\vec{x}\mathcal{L}_\vec{x}\phi = (-\mathcal{L}_\vec{y})(-\mathcal{L}_\vec{y})\phi = \mathcal{H}_\vec{y}\phi$.

$\blacksquare$

Construct the Interpolant

We wish to form an interpolant of the form $$ s(\vec{x}) = \sum c_i \phi(||\vec{x} - \vec{x}_i||) + d \mathcal{L}_\vec{y} \phi(||\vec{x} - \vec{y}||)+ e \mathcal{H}_\vec{y} \phi(||\vec{x} - \vec{y}||). $$ From Lemmas 1 and 3 we have that $\mathcal{L}_\vec{y} \phi(r) = -p \phi_1(r)$ and $\mathcal{H}_\vec{y} = \phi_1(r) + p^2 \phi_2(r)$.

We also require that at the stencil center $\vec{y}$ we have $\mathcal{H}_\vec{y}s(\vec{x}) \vert_{\vec{x}=\vec{y}} = 0$ and $\mathcal{L}_\vec{y}s(\vec{x}) \vert_{\vec{x}=\vec{y}} = 0$. Enforcing this will require us to evaluate the following six expressions: \begin{align*} &\mathcal{L}_\vec{x} \phi(||\vec{x}-\vec{x}_1||)|_{\vec{x}=\vec{y}} && \mathcal{L}_\vec{x}\mathcal{L}_\vec{y}\phi(r)|_{r=0} && \mathcal{L}_\vec{x}\mathcal{H}_\vec{y}\phi(r)|_{r=0} \\ &\mathcal{H}_\vec{x} \phi(||\vec{x}-\vec{x}_1||)|_{\vec{x}=\vec{y}} && \mathcal{H}_\vec{x}\mathcal{L}_\vec{y}\phi(r)|_{r=0} && \mathcal{H}_\vec{x}\mathcal{H}_\vec{y}\phi(r)|_{r=0}. \end{align*} From above we have that \begin{align*} \mathcal{L}_\vec{x}\mathcal{L}_\vec{y}\phi(r) &= -\mathcal{H}_\vec{x}\phi(r) \\ \mathcal{L}_\vec{x}\mathcal{H}_\vec{y}\phi(r) &= -\mathcal{L}_\vec{x}\mathcal{H}_\vec{x}\phi(r) \\ \mathcal{H}_\vec{x}\mathcal{L}_\vec{y}\phi(r) &= \mathcal{H}_\vec{x}\mathcal{L}_\vec{x}\phi(r) = \mathcal{L}_\vec{x}^3 = \mathcal{L}_\vec{x}\mathcal{H}_\vec{x}\phi(r) \\ \mathcal{H}_\vec{x}\mathcal{H}_\vec{y}\phi(r) &= \mathcal{H}_\vec{x}\mathcal{H}_\vec{x}\phi(r). \end{align*} and thus we need only calculate two expressions: $$ \mathcal{L}_\vec{x}\mathcal{H}_\vec{x}\phi(r) \phantom{=} \text{ and } \phantom{=} \mathcal{H}_\vec{x}\mathcal{H}_\vec{x}\phi(r). $$

We next observe that \begin{align*} \mathcal{L}_\vec{x}\mathcal{H}_\vec{x}\phi(r) &= \mathcal{L}_\vec{x} \big( \phi_1(r) + p^2 \phi_2(r) \big) \\ &= \mathcal{L}_\vec{x} \phi_1(r) + \phi_2(r)\mathcal{L}_\vec{x}p^2 + p^2 \mathcal{L}_\vec{x} \phi_2(r) \\ &= p\phi_2(r) + 2p \phi_2(r) + p^3 \phi_3(r) \\ &= 3p \phi_2(r) + p^3 \phi_3(r) \\ \mathcal{H}_\vec{x}\mathcal{H}_\vec{x} \phi(r) &= \mathcal{L}_\vec{x}\mathcal{L}_\vec{x}\mathcal{H}_\vec{x} \phi(r) \\ &= \mathcal{L}_\vec{x} \big( 3p\phi_2(r) + p^3 \phi_3(r) \big) \\ &= 3\phi_2(r)\mathcal{L}_\vec{x}p + 3p\mathcal{L}_\vec{x}\phi_2(r) + \phi_3(r) \mathcal{L}_\vec{x} p^3 + p^3 \mathcal{L}_\vec{x} \phi_3(r) \\ &= 3\phi_2(r) + 3p^2\phi_3(r) + 3p^2\phi_3(r) + p^4 \phi_4(r) \\ &= 3\phi_2(r) + 6p^2\phi_3(r) + p^4 \phi_4(r) \end{align*}

Finally we have our interpolation matrix equation $$ \begin{bmatrix} A & \vec{L} & \vec{H} \\ -\vec{L}^T & -\mathcal{H}_\vec{x}\phi(r)\vert_{r=0} & -\mathcal{L}_\vec{x}\mathcal{H}_\vec{x}\phi(r)\vert_{r=0} \\ \vec{H}^T & \mathcal{L}_\vec{x}\mathcal{H}_\vec{x}\phi(r)\vert_{r=0} & \mathcal{H}_\vec{x}\mathcal{H}_\vec{x} \phi(r)\vert_{r=0} \\ \end{bmatrix} \begin{bmatrix} \vec{c} \\ d \\ e \end{bmatrix} = \begin{bmatrix} \vec{f} \\ 0 \\0 \end{bmatrix} $$ where $A \in \RR^{n \times n}, \vec{L} \in \RR^n, \vec{H} \in \RR^n$ and \begin{align*} A_{ij} &= \phi(|| \vec{x}_i-\vec{x}_j ||) \\ \vec{L}_i &= \mathcal{L}_\vec{x}\phi(|| \vec{x}-\vec{x}_1 ||)\vert_{\vec{x} = \vec{x}_i} \\ \vec{H}_i &= \mathcal{H}_\vec{x}\phi(|| \vec{x}-\vec{x}_1 ||)\vert_{\vec{x} = \vec{x}_i}. \end{align*}

Take the Laplacian

Once the interpolant is formed it has the property that $\Delta s = \Delta_\SS s$ however we still must calculate $\Delta s$.To do so we must compute $\Delta \phi, \Delta\mathcal{L}_\vec{x}\phi$, and $\Delta\mathcal{H}_\vec{x}\phi$. From earlier work we have that $\Delta \phi(r) = 2\phi_1(r) + \phi''(r)$. In the notation of Lemma 1, let $\vec{x} = \begin{bmatrix} x & y & z \end{bmatrix}^T, \vec{y} = \begin{bmatrix} a & b & c \end{bmatrix}^T$, and $r = ||\vec{x} - \vec{y}||$. Additionally we will use the notation that Then we have that $p = \vec{n}\cdot(\vec{x}-\vec{y})$. \begin{align*} \frac{d}{dx}\mathcal{L}_\vec{x} \phi(r) &= \frac{d}{dx}p \phi_{1}(r) \\ &= n_x\phi_1(r) + p(x-a)\phi_2(r) \\ \frac{d^2}{dx^2}\mathcal{L}_\vec{x} \phi(r) &= n_x(x-a)\phi_2(r) + n_x(x-a)\phi_2(r) + p\phi_2(r) + p(x-a)^2 \phi_3(r) \\ &= 2n_x(x-a)\phi_2(r) + p\phi_2(r) + p(x-a)^2 \phi_3(r) \\ \Delta \mathcal{L}_\vec{x}\phi(r) &= 5p\phi_2(r) + r^2 p \phi_3(r). \end{align*} And lastly \begin{align*} \frac{d}{dx}\mathcal{H}_\vec{x} \phi(r) &= \frac{d}{dx} \big( \phi_1(r) + p^2 \phi_2(r) \big) \\ &= (x-a)\phi_2(r) + 2p n_x \phi_2(r) + p^2 (x-a) \phi_3(r) \\ \frac{d^2}{dx^2}\mathcal{H}_\vec{x} \phi(r) &= \phi_2(r) + (x-a)^2\phi_3(r) + 2n_x^2\phi_2(r) + 2p n_x(x-a) \phi_3(r) + \\ &\phantom{====} 2p n_x (x-a) \phi_3(r) + p^2 \phi_3(r) + p^2 (x-a)^2\phi_4(r) \\ &= \phi_2(r) + (x-a)^2\phi_3(r) + 2n_x^2\phi_2(r) + 4p n_x(x-a) \phi_3(r) + p^2 \phi_3(r) + p^2 (x-a)^2\phi_4(r) \\ \Delta \mathcal{H}_\vec{x} \phi(r) &= 3\phi_2(r) + r^2\phi_3(r) + 2\phi_2(r) + 4p^2 \phi_3(r) + 3p^2\phi_3(r) + p^2 r^2 \phi_4(r) \\ &= 5\phi_2(r) + \big( r^2 + 7p^2 \big) \phi_3(r) + p^2 r^2 \phi_4(r). \end{align*}

These values occur in our expression of $$ \Delta_\SS s(\vec{x}) = \Delta s(\vec{x}) = \sum c_i \Delta\phi(||\vec{x} - \vec{x}_i||) + d \Delta\mathcal{L}_x \phi(||\vec{x} - \vec{y}||)+ e \Delta\mathcal{H}_x \phi(||\vec{x} - \vec{y}||) $$ which in matrix form is $$ \begin{bmatrix} \Delta A & \Delta \vec{L} & \Delta \vec{H} \end{bmatrix} \begin{bmatrix} \vec{c} \\ d \\e \end{bmatrix} = \Delta \vec{f}. $$ Substituting our expression for the interpolant we have $$ \begin{bmatrix} \Delta A & \Delta \vec{L} & \Delta \vec{H} \end{bmatrix} \begin{bmatrix} A & \vec{L} & \vec{H} \\ -\vec{L}^T & -\mathcal{H}_\vec{x}\phi(r)\vert_{r=0} & -\mathcal{L}_\vec{x}\mathcal{H}_\vec{x}\phi(r)\vert_{r=0} \\ \vec{H}^T & \mathcal{L}_\vec{x}\mathcal{H}_\vec{x}\phi(r)\vert_{r=0} & \mathcal{H}_\vec{x}\mathcal{H}_\vec{x} \phi(r)\vert_{r=0} \\ \end{bmatrix}^{-1} \begin{bmatrix} \vec{f} \\ 0 \\0 \end{bmatrix} = \Delta \vec{f}. $$

References

1. Cécile Piret. The orthogonal gradients method: a radial basis functions method for solving partial differential equations on arbitrary surfaces. Journal of Computational Physics, 231(14):4662 – 4675, 2012. URL: http://www.sciencedirect.com/science/article/pii/S0021999112001477, doi:https://doi.org/10.1016/j.jcp.2012.03.007.


2. Jian-Jun Xu and Hong-Kai Zhao. An eulerian formulation for solving partial differential equations along a moving interface. Journal of Scientific Computing, 19(1):573–594, Dec 2003. URL: https://doi.org/10.1023/A:1025336916176, doi:10.1023/A:1025336916176.