3. 4. 2 Topological Neighborhood Considerations

The first step of a finite difference simulation is the determination of the neighborhood. In a very simple two-dimensional case, the neighborhood of a vertex is defined by the four vertices which are covered by the edges incident with the original vertex. It is also assumed that the respective vertex has to fulfill some geometric criterion to be considered as neighboring. In order to obtain all neighboring vertices of a given basis vertex, the traversal function $ VNV$ is introduced. Figures 3.11 and 3.12 show possible neighborhoods used for the finite difference simulation applied to structured grids and to unstructured meshes.

Figure 3.11: The finite difference scheme applied to a structured cell complex. Two neighborhoods $ \mathbf{b}$ are shown for the initial vertex $ \mathbf{a}$ . The neighborhood comprising five points can be used in order to determine differential operators up to second order. The neighborhood comprising nine points can be used in order to determine differential operators up to fourth order.
\includegraphics[width=10cm]{DRAWINGS/fvm_equi.eps}

Figure 3.12: The finite difference scheme applied to an unstructured cell complex. The neighborhood $ \mathbf{b}$ comprises five vertices $ \mathbf{c}$ . The vertices $ \mathbf{c}$ which are used for the evaluation of the respective differential operator can be either determined by finding the vertices with the smallest distance to the initial vertex $ \mathbf{a}$ or by determining the vertices which are incident with edges incident with the initial vertex $ \mathbf{a}$ .
\includegraphics[width=10cm]{DRAWINGS/fdm_topo.eps}

For each neighboring vertex $ \mathbf{w}$ of a basis vertex $ \mathbf{v}$ (also the vertex $ \mathbf{v}$ itself, which is within its own neighborhood), the following series expansion is written with respect to the location of the initial vertex $ \mathbf{v}$ :

$\displaystyle u(\mathbf{v}, \mathbf{w}) = u(\mathbf{v}) + \partial_x (x(\mathbf...
...f{v})) + \partial_y (y(\mathbf{w}) - x(\mathbf{v})) + \ldots + \mathcal{O}(h^n)$ (3.52)

For a given truncation error, an appropriate number of neighboring vertices has to be chosen. If a given truncation error is required for a Taylor expansion, for instance $ \mathcal{O}(h^3)$ , first order and second order terms have to be considered. Each partial derivative of $ u$ as well as $ u$ itself is considered as solution variable in a coordinate system. If a trunction error of higher order is required, the number of solution variables increases. Accordingly, a larger number of vertices (neighborhood) is required to solve the equation and determine the dependence of the position of the vertices on the value of the derivatives determined.

Moreover, it is possible to determine higher-order derivatives of a function by using a higher-order Taylor expansion. This implies that the use of higher-order derivatives necessarily increases the number of vertices (neighborhood) used for the Taylor expansion.

One obtains a local equation system which can be solved analytically in special cases, especially under the assumption that the local grid intervals are equal (which shall be assumed in the following). One obtains the following equation system

$\displaystyle \mathbf{u}(\mathbf{w})$ $\displaystyle = \mathcal{G}(\mathbf{w}) \cdot [\partial](\mathbf{w}) (+ [\mathcal{O}(h^N)])$    
$\displaystyle [\partial](\mathbf{w})$ $\displaystyle = \mathcal{G}(\mathbf{w})^{-1} \cdot \mathbf{u}(\mathbf{w}) (+ [\mathcal{O}(h^N)]) \; ,$ (3.53)

where the vector of function values $ u(\mathbf{w})$ is denoted as $ \mathbf{u}(\mathbf{w})$ . The vector $ [u, \partial_x u, \partial_y u, \ldots]$ is denoted as $ [\partial]$ evaluated on the vertex $ (\mathbf{w})$ . The matrix containing the geometrical coefficients is referred to as $ \mathcal{G}$ . By inverting the geometrical coefficient matrix $ \mathcal{G}$ , one obtains the vector of derivatives. The matrix $ \mathcal{G}$ can be written as follows:

$\displaystyle \mathcal{G} = \bigsqcup_{VNV}^{\underline{v}} [1; x - x_{\underli...
...erline{v}})(y - y_{\underline{v}}); (x - x_{\underline{v}})^2/2; \ldots]^T \; .$ (3.54)

The vector $ \mathbf{u}$ can be written as the vector of quantities within the neighborhood

$\displaystyle \mathbf{u} = \bigsqcup_{VNV} q = (q_1, \ldots, q_n)^T\; .$ (3.55)

A linear differential operator $ R = \mathcal{L}$ can be written as inner product of a given vector $ \mathbf{d}$ and the vector of derivatives $ [\partial]$ .

$\displaystyle R = \mathcal{L} = \langle \mathbf{d}, [\partial] \rangle ( = 0 )$ (3.56)

By extending the term of canonic partial differential operators $ [\partial]$ , one obtains the residual expression $ R$ . It has to be assured that the order of elements that are passed to the $ \bigsqcup_{VNV}$ operation is the same for the evaluation of the matrix $ \mathcal{G}$ and the vector of function values $ \mathbf{u}$ .

$\displaystyle R = \mathbf{d} \cdot \mathcal{G}^{-1} \cdot \mathbf{u}$ (3.57)

Inserting (3.54) without the second order terms and (3.55) into (3.56) yields the following expression.

$\displaystyle R = \mathbf{d} \cdot [\bigsqcup_{VNV}^{\underline{v}} [1; x - x_{\underline{v}}; y - y_{\underline{v}}; \ldots]^T]^{-1} \cdot \bigsqcup_{VNV} q$ (3.58)

If the geometric constants are known explicitly, the inverse of the matrix $ \mathcal{G}$ can be derived directly. In this case the matrix $ \mathcal{G}$ can be written as:

$\displaystyle \mathcal{G} =
\left[\begin{array}{c c c c c}
1 & 0 & 0 & 0 & 0 \\...
...2 & 0 \\
1 & 0 & h & 0 & h^2/2 \\
1 & 0 & -h & 0 & h^2/2
\end{array}\right]
$

The inverse $ \mathcal{G}^{-1}$ of the matrix $ \mathcal{G}$ can be obtained as follows:

$\displaystyle \mathcal{G}^{-1} =
\left[\begin{array}{c c c c c}
1 & 0 & 0 & 0 &...
...h^2 & 1/h^2 & 0 & h^2/2 \\
-2/h^2 & 0 & 0 & 1/h^2 & 1/h^2
\end{array}\right]
$

The standard finite difference formulae can be obtained easily from the coefficients of the matrix. By the specification of $ \mathbf{d}$ , a linear combination of lines of the matrix $ \mathcal{G}^{-1}$ is obtained by multiplication. The line vector $ \mathbf{k} = \mathbf{d} \cdot \mathcal{G}^{-1}$ denotes the coefficients with which the values of the single neighboring vertices are coupling. It is associated with a vertex and contains coupling coefficients, each of which is directly associated with a neighboring vertex. Alternatively, the vector $ \mathbf{k}(\mathbf{w})$ can be written element-wise, where the elements are written as $ K(\mathbf{w}, \mathbf{v})$ . In this case $ \mathbf{w}$ denotes the vertex and $ \mathbf{v}$ denotes a vertex in the neighborhood of $ \mathbf{w}$ . By evaluating the inner product of $ \mathbf{k}$ with the function values of the single vertices stored in the vector $ \mathbf{u}$ , the residual expression is obtained. The evaluation of the inner product finally yields

$\displaystyle R = [\mathbf{k} \cdot \bigsqcup_{VNV} q] = [\bigsqcup_{VNV}^{\und...
...t \bigsqcup_{VNV} q] = [\sum_{VNV}^{\underline{w}} K(\underline{w}, \bullet) q]$ (3.59)

It can be seen easily that the algebraic structure of (3.59) is similar to other discretization schemes. While other discretization schemes use topologically specified double sums with topological traversal in order to obtain the couplings between solution quantities, the definition of the neighborhood can be defined more freely. While in one case the neighborhood is defined topologically like for finite elements or finite volumes, the neighborhood can also be determined by geometrical considerations.

Michael 2008-01-16