<figure><img src="image_1.png" /></figure>

HAL
open science

# Green's functions, boundary elements and finite elements

Friedel Hartmann

► To cite this version:

Friedel Hartmann. Green's functions, boundary elements and finite elements. 11th International Conference on Advances in Boundary Element Techniques, Jul 2010, Berlin, Germany. hal-01586562

HAL Id: hal-01586562

https://hal.science/hal-01586562

Submitted on 13 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Green's functions, boundary elements and finite elements

Friedel Hartmann

University of Kassel, Kurt-Wolters-Str. 3, D-34225

friedelhartmann@uni-kassel.de

Keywords: Green's function, boundary elements, finite elements, stiffness matrices, adaptive refinement

## Introduction

The Green's function $G(y, x)$ of the Poisson equation

$$ -\Delta u = p \quad \text{in } \Omega \quad u = 0 \quad \text{on } \Gamma \qquad (1) $$

allows to write the solution $u(x)$ in terms of an influence function

$$ u(x) = \int_{\Omega} G(y, x) p(y) \, d\Omega y. \qquad (2) $$

In boundary element analysis we replace the Green's function by a fundamental solution $g(y, x)$

$$ u(x) = \int_{\Gamma} [g(y, x) \frac{\partial u}{\partial n}(y) - \frac{\partial g}{\partial n}(y, x) u(y)] \, dsy + \int_{\Omega} g(y, x) p(y) \, d\Omega y \qquad (3) $$

and we determine the unknown boundary values by solving an integral equation approximately and so construct the BE-solution

$$ u_h(x) = \int_{\Gamma} [g(y, x) \frac{\partial u_h}{\partial n}(y) - \frac{\partial g}{\partial n}(y, x) u_h(y)] \, dsy + \int_{\Omega} g(y, x) p(y) \, d\Omega y. \qquad (4) $$

Solving the problem (1) with finite elements means that we project the exact solution $u(x)$ onto the solution and trial space

$$ V_h = \{v_h \in H_1(\Omega) | v_h(x) = 0, x \in \Gamma\} \qquad (5) $$

that is the FE-solution $u_h$ is the best approximation in the sense of the strain energy metric

$$ a(u - u_h, u - u_h) \leq a(u - v_h, u - v_h) \quad v_h \in V_h \qquad (6) $$

where

$$ a(u, v) := \int_{\Omega} (\nabla u \cdot \nabla v) \, d\Omega. \qquad (7) $$

## Tottenham's equation

It seems that these two methods follow two totally different lines of reasoning but the surprising result is that also the FE-method is a Green's function method. This result goes back to Tottenham who in a very early paper, [1], mentioned that the FE-solution $u_h$ has the form

$$ u_h(x) = \int_{\Omega} G_h(y, x) p(y) \, d\Omega y \qquad (8) $$

where $G_h(y, x)$ is the FE-Green's function that is the projection of the exact Green's function

$$ -\Delta y G(y, x) = \delta(y - x) \quad G(y, x) = 0 \quad y \in \Gamma \qquad (9) $$
onto the space $V_h$. If the nodal shape functions $\varphi_i(\mathbf{x})$ form a basis of $V_h$ then $G_h(\mathbf{y}, \mathbf{x})$ is the function

$$G_h(\mathbf{y}, \mathbf{x}) = \sum_i u_i(\mathbf{x}) \varphi_i(\mathbf{y}) \qquad (10)$$

where the nodal values $u_i = u_i(\mathbf{x})$ are the solution of the system

$$K \mathbf{u}(\mathbf{x}) = \mathbf{f} \quad k_{ij} = a(\varphi_i, \varphi_j) \quad f_i = \varphi_i(\mathbf{x}). \qquad (11)$$

The proof of Tottenham's equation (8) is easy and can be done on one line. Namely we have

$$u_h(\mathbf{x}) = \int_{\Omega} \delta(\mathbf{y} - \mathbf{x}) u_h(\mathbf{y}) d\Omega \mathbf{y} = a(G_h, u_h) = \int_{\Omega} G_h(\mathbf{y}, \mathbf{x}) p(\mathbf{y}) d\Omega \mathbf{y}. \qquad (12)$$

First $u_h$ is considered to be the FE solution of the right-hand side $p$ and $G_h \in V_h$ assumes the role of a virtual displacement

$$a(G_h, u_h) = \int_{\Omega} G_h(\mathbf{y}, \mathbf{x}) p(\mathbf{y}) d\Omega \mathbf{y} \qquad (13)$$

next $G_h$ is considered to be the FE solution of the right-hand side $\delta$, and $u_h$ assumes the role of a virtual displacement

$$u_h(\mathbf{x}) = \int_{\Omega} \delta(\mathbf{y} - \mathbf{x}) u_h(\mathbf{y}) d\Omega \mathbf{y} = a(G_h, u_h) \qquad (14)$$

which explains the left-hand side. The symmetric strain energy $a(G_h, u_h)$ plays the role of a turnstile.

## The inverse stiffness matrix

The approximate Green's function for the displacement $u(\mathbf{x})$ at a node $\mathbf{x}_k$ has the form

$$G_h(\mathbf{y}, \mathbf{x}_k) = \sum_i u_{Gi}(\mathbf{x}_k) \varphi_i(\mathbf{y}). \qquad (15)$$

The vector $\mathbf{u}_G = \{u_{G1}, u_{G2}, \dots, u_{Gn}\}$ is the solution of the $n \times n$ system

$$K \mathbf{u}_G = \mathbf{e}_k \quad (\text{unit vector } \mathbf{e}_k), \qquad (16)$$

which means that the columns $\mathbf{c}_k$ of the inverse stiffness matrix $K^{-1}$

$$\mathbf{u}_G = K^{-1} \mathbf{e}_k = \mathbf{c}_k \qquad (17)$$

are the nodal displacements which belong to the $n$ Green's functions $G_h(\mathbf{y}, \mathbf{x}_k)$ of the $n$ nodes $\mathbf{x}_k$.

$$G_h(\mathbf{y}, \mathbf{x}_k) = \sum_i c_{ki} \varphi_i(\mathbf{y}) = \mathbf{c}_k^T \varphi_i(\mathbf{y}). \qquad (18)$$

This explains why the inverse of a tridiagonal matrix is fully populated. Even if only one node $\mathbf{x}_k$ carries a point load $P = 1$ the whole structure deforms. An example provides the bar in Fig. 1a, which consists of five linear elements. The stiffness matrix and its inverse are, see Fig. 1b

$$K = \frac{EA}{l} \begin{bmatrix} 2 & -1 & 0 & 0 \\ -1 & 2 & -1 & 0 \\ 0 & -1 & 2 & -1 \\ 0 & 0 & -1 & 2 \end{bmatrix}, \quad \Rightarrow \quad K^{-1} = \frac{l}{5EA} \begin{bmatrix} 4 & 3 & 2 & 1 \\ 3 & 6 & 4 & 2 \\ 2 & 4 & 6 & 3 \\ 1 & 2 & 3 & 4 \end{bmatrix}. \qquad (19)$$

And evidently are the columns of the inverse the Green's functions of the nodes.
<figure><img src="image_3.png" /><figcaption>Figure 1: a) Elastic bar subdivided into five linear elements b-e) the displacements are the columns of the inverse stiffness matrix (all values times l/(5 EA)).</figcaption></figure>

## Boundary elements

When the BE solution (4) is compared with the FE solution (8) there seems not to be much agreement between the two solutions. But the BE solution (4) is just the same formula

$$u_h(\mathbf{x}) = \int_{\Omega} G_h(\mathbf{y}, \mathbf{x}) p(\mathbf{y}) \, d\Omega \mathbf{y} \qquad (20)$$

in disguise, of course with a different $G_h$. To see this note that the Green's function G can be split into a fundamental solution G and a regular part $u_R$,

$$u(\mathbf{x}) = \int_{\Omega} G(\mathbf{y}, \mathbf{x}) p(\mathbf{y}) \, d\Omega \mathbf{y} = \int_{\Omega} g(\mathbf{y}, \mathbf{x}) p(\mathbf{y}) \, d\Omega \mathbf{y} + \int_{\Omega} u_R(\mathbf{y}, \mathbf{x}) p(\mathbf{y}) \, d\Omega \mathbf{y}. \qquad (21)$$

The regular part is the solution of the problem

$$-\Delta u_R = 0 \quad \text{in } \Omega \quad u_R(\mathbf{y}, \mathbf{x}) = -g(\mathbf{y}, \mathbf{x}) \quad \mathbf{y} \in \Gamma. \qquad (22)$$

This splitting implies that the boundary integrals in the influence function (4) are just an equivalent expression for the work done by the distributed load p acting through the regular part $u_R$:

$$\int_{\Omega} u_R(\mathbf{y}, \mathbf{x}) p(\mathbf{y}) \, d\Omega \mathbf{y} = \int_{\Gamma} G(\mathbf{y}, \mathbf{x}) \frac{\partial u}{\partial n}(\mathbf{y}) \, ds \mathbf{y} - \int_{\Gamma} \frac{\partial G}{\partial \nu}(\mathbf{y}, \mathbf{x}) u(\mathbf{y}) \, ds \mathbf{y}. \qquad (23)$$

Hence it can be assumed that the boundary integrals in the BE solution (4) play the same role,

$$\int_{\Omega} u_R^h(\mathbf{y}, \mathbf{x}) p(\mathbf{y}) \, d\Omega \mathbf{y} := \int_{\Gamma} G(\mathbf{y}, \mathbf{x}) \frac{\partial u_h}{\partial n}(\mathbf{y}) \, ds \mathbf{y} - \int_{\Gamma} \frac{\partial G}{\partial \nu}(\mathbf{y}, \mathbf{x}) u_h(\mathbf{y}) \, ds \mathbf{y}, \qquad (24)$$
and so we arrive at (20), which makes the two methods look alike. Though there is a difference in the Green's function: the FE Green's function is only mesh dependent, while the BE Green's function is also load-case dependent.

Note that the BE method uses an ingenious approach: it does not approximate $u_R(y, x)$—this would be simply too laborious, because at every point $x$ it would have to approximate a new function $u_R(y, x)$—but instead substitutes for the work integral $(u_R[x], p)$ the work done by the Cauchy data $u$ on $\Gamma_N$ and $t = \partial u/\partial n$ on $\Gamma_D$ via the conjugate terms of the fundamental solution $G$. That is, the program knows that it suffices to approximate the “static” data $u$ and $t$ leaving the effects caused by a change in the observation point $x$ to the fundamental solution $G(y, x)$. This is the essence of (24).

## Accuracy

The FE method has a handicap because it must approximate both parts of the Green's function $G = g + u_R$ while in the BE method the fundamental solution is built into the code. This handicap can be overcome if the FE-method approximates only the regular part

$$ 
\begin{aligned}
u(\boldsymbol{x}) - u_h(\boldsymbol{x}) &= \int_{\Omega} (G(y, \boldsymbol{x}) - G_h(y, \boldsymbol{x})) p(\boldsymbol{x}) d\Omega y \\
&= \int_{\Omega} (g(y, \boldsymbol{x}) + u_R(\boldsymbol{x}) - g(y, \boldsymbol{x}) + u_R^h(y, \boldsymbol{x})) p(y) d\Omega y \\
&= \int_{\Omega} (u_R(\boldsymbol{x}) - u_R^h(y, \boldsymbol{x})) p(y) d\Omega y.
\end{aligned} 
\quad (25) $$

Of course if the source point $x$ comes too close to the boundary then the regular part will also become singular.

## Everything is a functional

The importance of Tottenham's insight lies in the fact that in linear analysis every value, every displacement, every stress, every potential value can be considered a functional $J(u)$. For example the following expressions

$$ J(u) = u(\boldsymbol{x}) \qquad J(u) = \sigma_{xx}(\boldsymbol{x}) \qquad J(u) = \int_{0}^{l} \sigma_{xy} ds \quad (26) $$

are all functionals. And so as we associate with the classical functional, $J(u) = u(\boldsymbol{x})$, a Dirac delta

$$ J(u) = \int_{\Omega} \delta(y - \boldsymbol{x}) u(y) d\Omega y \quad (27) $$

so we can associate with every other linear functional $J(u)$ a certain Dirac delta. Let the functional $J(u)$ be for example the average value of the gradient of $u$ over a small circular region $\Omega_c(\boldsymbol{x})$ with radius $c$ centered at a point $\boldsymbol{x} \in \Omega$

$$ J(u) = \frac{1}{\pi c^2} \int_{\Omega_c(\boldsymbol{x})} (u_{,x}(\boldsymbol{x}) + u_{,y}(\boldsymbol{x})) d\Omega. \quad (28) $$

We postulate that this functional can be expressed as

$$ J(u) = \int_{\Omega} \delta_{\Sigma}(y - \boldsymbol{x}) u(y) d\Omega y. \quad (29) $$

We do not say what the Dirac delta $\delta_{\Sigma}(y - \boldsymbol{x})$ looks like, what it is. We only say that the scalar product (= integral) of $\delta_{\Sigma}$ with the function $u$ provides the value $J(u)$. In more complicated cases it
is even questionable whether the Dirac delta is an integrable function, whether the expressions make sense. We simply assume they do.

Next we only need to find the function G which has the delta function as its right-hand side, that is which solves the equation

$$ -\Delta G(\mathbf{y}, \mathbf{x}) = \delta_{\Sigma}(\mathbf{y} - \mathbf{x}) \quad G = 0 \quad \text{on } \Gamma \qquad (30) $$

because Green's second identity then implies

$$ B(G, u) = \underbrace{\int_{\Omega} \delta_{\Sigma}(\mathbf{y} - \mathbf{x}) u(\mathbf{y}) d\Omega}_{J(u)} - \int_{\Omega} G(\mathbf{x}, \mathbf{y}) p(\mathbf{y}) d\Omega_{\mathbf{y}} = 0 \qquad (31) $$

or

$$ J(u) = \int_{\Omega} G(\mathbf{x}, \mathbf{y}) p(\mathbf{y}) d\Omega_{\mathbf{y}}. \qquad (32) $$

This is the 'trick'. The abstract Dirac delta has got a shape, G, and the scalar product of this G with the right hand side p renders the value J(u).

Now we do finite elements. We approximate the Green's function G with the n shape functions in $V_h$

$$ G_h(\mathbf{y}, \mathbf{x}) = \sum_{i=1}^{n} u_i(\mathbf{x}) \varphi_i(\mathbf{y}) \qquad (33) $$

and so the variational problem for the weak solution becomes on $V_h$ a system of n equations

$$ a(G_h, \varphi_i) = (\delta_{\Sigma}(\mathbf{y} - \mathbf{x}), \varphi_i) = J(\varphi_i) \quad i = 1, 2, \dots, n \qquad (34) $$

or in matrix notation

$$ K u(\mathbf{x}) = f(\mathbf{x}) \quad k_{ij} = a(\varphi_i, \varphi_j) \quad f_i = J(\varphi_i). \qquad (35) $$

So the extension of the FE method to Green's functions is simple and straightforward.

## Infinite energy

Strictly speaking the extension is not so simple because most Green's functions have infinite energy. Hence approximating a Green's function with finite elements is an ill-posed problem. The surprising aspect is that the FE-method does not care. All the input in an FE-program is processed by these approximate Green's function which theoretically do not lie in the solution space.

But it is remarkable how 'old-fashioned' the ansatz of the FE-method is. The FE-method is a numerical implementation of Green's function. And this is the same idea behind the boundary element method. So in this sense the two methods are identical twins.

## Goal oriented refinement

The error of an FE-solution

$$ u(\mathbf{x}) - u_h(\mathbf{x}) = \int_{\Omega} (G(\mathbf{y}, \mathbf{x}) - G_h(\mathbf{y}, \mathbf{x})) p(\mathbf{x}) d\Omega_{\mathbf{y}} \qquad (36) $$

is simply proportional to the error in the approximate Green's function. So given that we are interested in the value $u(\mathbf{x})$ of the solution at a particular point $\mathbf{x}$ we must refine the mesh in such a way that the error in the Green's function for the point value $u(\mathbf{x})$ is minimized. This technique is called *goal oriented adaptive refinement*.
By making use of the Galerkin orthogonality and Schwarz' inequality we can formulate the following estimate for the error of an FE-solution

$$|u(\mathbf{x}) - u_h(\mathbf{x})| \le ||\mathbf{G}[\mathbf{x}] - \mathbf{G}_h[\mathbf{x}]||_E ||\mathbf{u} - \mathbf{u}_h||_E \quad (37)$$

where $||\cdot||_E$ is the standard energy norm, $||\mathbf{u}||_E^2 := a(\mathbf{u}, \mathbf{u})$ and $G[\mathbf{x}] := G(\mathbf{y}, \mathbf{x})$ is the Green's function. Again we face the same problem as before: theoretically such an estimate requires the Green's function to have bounded energy what normally is not the case. But if we read $u(\mathbf{x})$ as the average value of the potential over a small disk centered at $\mathbf{x}$ the corresponding Green's function $\mathbf{G}(\mathbf{y}, \mathbf{x})$ would have finite energy and then (37) would make sense.

This result (37) is motivation to minimize the error in the Green's function and the error in the solution $u_h$ simultaneously. That is, at each refinement step an error indicator $\eta_G^e$ for the Green's function and an error indicator $\eta_p^e$ for the original problem is calculated on each element and the combined error indicator on each element is then $\eta^e = \eta_G^e \cdot \eta_p^e$ so that the sum of the local errors $\eta^e = \eta_G^e \cdot \eta_p^e$ provides an upper bound for the error

$$|u(\mathbf{x}) - u_h(\mathbf{x})| \le \sum_e \eta_G^e \cdot \eta_p^e. \quad (38)$$

The energy norms of the errors $e_G = G - G_h$ and $e_u = u - u_h$ are calculated by measuring the eigenwork done by the element residual forces and jump terms on the element edges, so that for example

$$\eta_p = a(e_u, e_u)_{\Omega_e}^{1/2} = r(e_u)^{1/2}. \quad (39)$$

This technique was applied to the plate in Fig. 2. The first mesh, the mesh in Fig. 2 a, is the result of a standard adaptive refinement

$$\eta_p^e \le \varepsilon_{TOL}. \quad (40)$$

To push the error below the preset error margin the program has to refine all those elements—in practice only the first, say, 30%—where the error $\eta_p^e$ of the original (or primal) problem exceeds this margin.

The mesh in Fig. 2 b is based on weighting the primal error $\eta_p$ with the error $\eta_G$ of the Green's function elementwise

$$\eta_G^e \cdot \eta_p^e \le \varepsilon_{TOL}. \quad (41)$$

In most parts of the mesh the error $\eta_G^e$ of the numerical Green's function is so low that this inequality is automatically satisfied. That is many of the refinements in Fig. 2 a are not necessary if we are only interested in the stress $\sigma_{yy}$ at $\mathbf{x}$.

## References

[1] Tottenham H (1970) "Basic Principles", in: Finite Element Techniques in Structural Mechanics. (Eds. Tottenham H, Brebbia C), Southampton University Press, Southampton 1970

[2] Hartman F, Katz C (2008) Structural Analysis with Finite Elements, 2nd ed., Springer-Verlag Berlin Heidelberg New York
<figure><img src="image_4.png" /><figcaption>Figure 2: Adaptive refinement **a)** standard refinement $\eta_p \le \varepsilon_{TOL}$ **b)** goal oriented refinement $\eta_p \times \eta_G \le \varepsilon_{TOL}$</figcaption></figure>