9.2.5. Transformation auf das Referenzelement

In der Finite-Elemente-Methode werden Basisfunktionen und Quadraturregeln nur einmal auf einem Referenzelement \(T\) definiert. Jedes Element \(T_i\) der Triangulierung wird über eine affine Abbildung \(\sigma_i\) darauf zurückgeführt.

Allgemeines Prinzip

Gegeben sei ein nicht entartetes Simplex \(T_i\subset\mathbb{R}^d\) mit Ecken \(\vec{p}_{i,1},\dots,\vec{p}_{i,d+1}\) sowie das Referenzsimplex

\[ T = \Bigl\{\vec{t}\in\mathbb{R}^d \;:\; t_k\ge 0,\ \textstyle\sum_{k=1}^d t_k \le 1\Bigr\} \]

mit den Ecken \(\vec{0}\) und \(\vec{e}_k\) (\(k=1,\dots,d\)). Die Koordinaten im Referenzelement heissen \(\xi\) (1D), \((\xi,\eta)\) (2D) bzw. \((\xi,\eta,\zeta)\) (3D), zusammengefasst im Vektor \(\vec{t}\).

Abbildung Referenz → Element. Die affine Abbildung

(9.14)\[\sigma_i : T \to T_i,\qquad \vec{x} = \sigma_i(\vec{t}\,) = \vec{p}_{i,1} + A_i\cdot\vec{t}, \qquad A_i = \bigl(\vec{p}_{i,2}-\vec{p}_{i,1},\ \vec{p}_{i,3}-\vec{p}_{i,1},\ \dots,\ \vec{p}_{i,d+1}-\vec{p}_{i,1}\bigr)\in\mathbb{R}^{d\times d}\]

bildet die Ecken des Referenzelements auf die Ecken \(\vec{p}_{i,k}\) ab. Die Spalten der Jacobi-Matrix \(A_i = D\sigma_i\) sind die von \(\vec{p}_{i,1}\) ausgehenden Kantenvektoren. Die Jacobi-Determinante bezeichnen wir mit

\[ J_i = \det A_i . \]

Rücktransformation Element → Referenz. Ist \(J_i\neq 0\) (d.h. \(T_i\) ist nicht entartet), so gilt

(9.15)\[\sigma_i^{-1} : T_i \to T,\qquad \vec{t} = \sigma_i^{-1}(\vec{x}) = A_i^{-1}\cdot\bigl(\vec{x}-\vec{p}_{i,1}\bigr), \qquad D\sigma_i^{-1}(\vec{x}) = A_i^{-1}.\]

Wichtige Folgerungen.

• Volumen: \(|T_i| = |J_i|\,|T|\) mit \(|T| = 1/d!\).

• Integrale: \(\displaystyle\int_{T_i} f(\vec{x})\,d\vec{x} = \int_{T} f\bigl(\sigma_i(\vec{t}\,)\bigr)\,|J_i|\,d\vec{t}\ \).

• Gradienten: Für \(\varphi(\vec{x})=\tilde{\varphi}(\sigma_i^{-1}(\vec{x}))\) gilt \(D\varphi(\vec{x}) = D\tilde{\varphi}(\sigma_i^{-1}(\vec{x}))\cdot A_i^{-1}\), also \(\nabla\varphi(\vec{x}) = {A_i^{-1}}^T\cdot\nabla\tilde{\varphi}(\sigma_i^{-1}(\vec{x}))\).

• Orientierung: \(J_i>0\) genau dann, wenn \(\vec{p}_{i,1},\dots,\vec{p}_{i,d+1}\) gleich orientiert sind wie die Referenzecken.

• Die Referenzkoordinaten \(\vec{t}\) sind eng mit den baryzentrischen Koordinaten von \(\vec{x}\) verwandt, siehe Baryzentrische Koordinaten.

Baryzentrische Koordinaten

Definition. Sei \(T_i\) ein nicht entartetes Simplex mit Ecken \(\vec{p}_{i,1},\dots,\vec{p}_{i,d+1}\). Jeder Punkt \(\vec{x}\in\mathbb{R}^d\) lässt sich eindeutig als gewichtete Summe der Ecken schreiben, deren Gewichte sich zu eins summieren:

(9.16)\[\vec{x} = \sum_{k=1}^{d+1}\lambda_k\,\vec{p}_{i,k},\qquad \sum_{k=1}^{d+1}\lambda_k = 1 .\]

Die Zahlen \(\lambda_1(\vec{x}),\dots,\lambda_{d+1}(\vec{x})\) heissen baryzentrische Koordinaten von \(\vec{x}\) bezüglich \(T_i\).

Zusammenhang mit den Referenzkoordinaten. Setzt man \(\lambda_1 = 1-\sum_{k=2}^{d+1}\lambda_k\) in (9.16) ein, so erhält man

\[\begin{split} \vec{x}-\vec{p}_{i,1} = \sum_{k=2}^{d+1}\lambda_k\,\bigl(\vec{p}_{i,k}-\vec{p}_{i,1}\bigr) = A_i\cdot\begin{pmatrix}\lambda_2\\\vdots\\\lambda_{d+1}\end{pmatrix}. \end{split}\]

Das ist genau die Gleichung \(\vec{x}=\sigma_i(\vec{t}\,)\). Weil \(A_i\) invertierbar ist, sind die \(\lambda_k\) eindeutig bestimmt, und es gilt

\[ \lambda_{k+1}(\vec{x}) = t_k \quad (k=1,\dots,d),\qquad \lambda_1(\vec{x}) = 1-\sum_{k=1}^{d} t_k, \qquad\text{mit}\quad \vec{t}=\sigma_i^{-1}(\vec{x}). \]

Im Dreieck also \(\lambda_1 = 1-\xi-\eta\), \(\lambda_2=\xi\), \(\lambda_3=\eta\).

Eigenschaften.

• An den Ecken gilt \(\lambda_k(\vec{p}_{i,j}) = \delta_{kj}\). Da die \(\lambda_k\) affin in \(\vec{x}\) sind, sind sie genau die linearen Basisfunktionen auf \(T_i\): \(\lambda_{k+1}(\vec{x}) = \tilde{\varphi}_k(\sigma_i^{-1}(\vec{x}))\) mit \(\tilde{\varphi}_0 = 1-\xi-\eta\), \(\tilde{\varphi}_1=\xi\), \(\tilde{\varphi}_2=\eta\) (2D).

• \(\vec{x}\in T_i\) genau dann, wenn alle \(\lambda_k(\vec{x})\ge 0\). Auf der Seite gegenüber von \(\vec{p}_{i,k}\) ist \(\lambda_k=0\).

• Der Schwerpunkt von \(T_i\) hat die baryzentrischen Koordinaten \(\lambda_k=\frac{1}{d+1}\).

Geometrische Deutung. Ersetzt man im Simplex \(T_i\) die Ecke \(\vec{p}_{i,k}\) durch den Punkt \(\vec{x}\), so entsteht ein Teilsimplex \(T_i^{(k)}(\vec{x})\). Es gilt

\[ \lambda_k(\vec{x}) = \frac{|T_i^{(k)}(\vec{x})|}{|T_i|}, \]

wobei die Volumina orientiert zu nehmen sind (für \(\vec{x}\notin T_i\) werden einzelne \(\lambda_k\) negativ). Das folgt aus der Cramerschen Regel für \(A_i\cdot\vec{t}=\vec{x}-\vec{p}_{i,1}\): Ersetzt man in \(A_i\) eine Spalte durch \(\vec{x}-\vec{p}_{i,1}\), so ist die Determinante dieser Matrix das \(d!\)-fache Volumen des entsprechenden Teilsimplex. Für \(\lambda_1\) folgt die Aussage dann aus \(\sum_k\lambda_k=1\), da sich die Teilvolumina zu \(|T_i|\) summieren. Im Dreieck ist \(\lambda_k\) also der Anteil des Teildreiecks, das der Ecke \(\vec{p}_{i,k}\) gegenüberliegt, an der Fläche von \(T_i\); im Tetraeder entsprechend der Volumenanteil eines Teiltetraeders.

Abb. 9.8 Baryzentrische Koordinaten im Dreieck: \(\lambda_k(\vec{x})\) ist der Flächenanteil des Teildreiecks, das der Ecke \(\vec{p}_{i,k}\) gegenüberliegt. Im Beispiel ist \(\vec{x} = 0.2\,\vec{p}_{i,1} + 0.5\,\vec{p}_{i,2} + 0.3\,\vec{p}_{i,3}\).

1D: Interval

Mit \(p_{i,1}=a\) und \(p_{i,2}=b\) ist \(A_i=J_i=b-a\) ein Skalar. Die Referenzkoordinate ist hier der Skalar \(\vec{t}=\xi\in[0,1]\).

\[\begin{split} \begin{aligned} \text{Abbildung:}\quad & x = \sigma_i(\xi) = a + (b-a)\,\xi, \\ \text{Rücktransformation:}\quad & \xi = \sigma_i^{-1}(x) = \frac{x-a}{b-a}, \\ \text{Jacobi:}\quad & J_i=\frac{dx}{d\xi} = b-a,\qquad \frac{d\xi}{dx}=\frac{1}{b-a},\\ \text{Integral:}\quad & \int_a^b f(x)\,dx = (b-a)\int_0^1 f\bigl(a+(b-a)\xi\bigr)\,d\xi,\\ \text{Ableitung:}\quad & \varphi'(x) = \frac{1}{b-a}\,\tilde{\varphi}'\bigl(\sigma_i^{-1}(x)\bigr). \end{aligned} \end{split}\]

Beispiel. \([a,b]=[1,4]\): \(J_i=3\), \(\xi=(x-1)/3\), \(x = 1+3\xi\). Der Mittelpunkt \(x=2.5\) wird auf \(\xi = 0.5\) abgebildet.

Abb. 9.9 Affine Abbildung zwischen \(T_i=[a,b]\) und dem Einheitsintervall \(T=[0,1]\).

2D: Dreieck

Seien \(\vec{p}_{i,k}=(x_k,y_k)^T\) die Ecken von \(T_i\). Das Referenzdreieck \(T\) hat die Ecken \((0,0)\), \((1,0)\) und \((0,1)\). Ein Punkt in \(T\) wird mit \(\vec{t}=(\xi,\eta)^T\) bezeichnet.

Abbildung.

\[\begin{split} \sigma_i(\vec{t}\,) = \sigma_i(\xi,\eta) = \vec{p}_{i,1} + A_i\cdot\begin{pmatrix}\xi\\\eta\end{pmatrix}, \qquad A_i = \begin{pmatrix} x_2-x_1 & x_3-x_1\\ y_2-y_1 & y_3-y_1\end{pmatrix}, \end{split}\]

bzw. in baryzentrischer Form \(\sigma_i(\xi,\eta) = (1-\xi-\eta)\,\vec{p}_{i,1} + \xi\,\vec{p}_{i,2} + \eta\,\vec{p}_{i,3}\).

Jacobi-Determinante.

\[ J_i = \det A_i = (x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1) = \pm 2\,|T_i| , \qquad dx\,dy = |J_i|\,d\xi\,d\eta . \]

\(J_i\) ist genau dann positiv, wenn \(\vec{p}_{i,1},\vec{p}_{i,2},\vec{p}_{i,3}\) im Gegenuhrzeigersinn angeordnet sind.

Inverse Jacobi-Matrix.

\[\begin{split} A_i^{-1} = \frac{1}{J_i} \begin{pmatrix} y_3-y_1 & -(x_3-x_1)\\ -(y_2-y_1) & x_2-x_1\end{pmatrix}. \end{split}\]

Rücktransformation.

\[\begin{split} \begin{pmatrix}\xi\\\eta\end{pmatrix} = \sigma_i^{-1}(x,y) = A_i^{-1}\cdot\left(\begin{pmatrix}x\\y\end{pmatrix}-\vec{p}_{i,1}\right), \qquad \begin{aligned} \xi &= \frac{(y_3-y_1)(x-x_1)-(x_3-x_1)(y-y_1)}{J_i},\\ \eta &= \frac{(x_2-x_1)(y-y_1)-(y_2-y_1)(x-x_1)}{J_i}. \end{aligned} \end{split}\]

Nach Baryzentrische Koordinaten ist \(\xi=\lambda_2(\vec{x})\) und \(\eta=\lambda_3(\vec{x})\): Der Zähler von \(\xi\) ist das Doppelte der orientierten Fläche des Teildreiecks \((\vec{p}_{i,1},\vec{x},\vec{p}_{i,3})\), der Zähler von \(\eta\) jenes von \((\vec{p}_{i,1},\vec{p}_{i,2},\vec{x})\). Geteilt durch \(J_i = 2\,|T_i|\) ergeben sich die Flächenanteile.

Beispiel. \(\vec{p}_{i,1}=(1,1)^T\), \(\vec{p}_{i,2}=(4,2)^T\), \(\vec{p}_{i,3}=(2,4)^T\):

\[\begin{split} A_i=\begin{pmatrix}3&1\\1&3\end{pmatrix},\quad J_i = 8,\quad |T_i| = 4,\quad A_i^{-1}=\frac18\begin{pmatrix}3&-1\\-1&3\end{pmatrix}. \end{split}\]

Probe: \(\sigma_i^{-1}(\vec{p}_{i,2})=\tfrac18(3\cdot 3-1,\,-3+3)^T=(1,0)^T\) und \(\sigma_i^{-1}(\vec{p}_{i,3})=(0,1)^T\). Der Schwerpunkt \((\tfrac73,\tfrac73)^T\) wird auf \((\tfrac13,\tfrac13)^T\) abgebildet.

Abb. 9.10 Affine Abbildung zwischen einem beliebigen Dreieck \(T_i\) und dem Einheitsdreieck \(T\). Die grauen Linien \(\xi=\text{const}\) und \(\eta=\text{const}\) werden auf parallele Geraden abgebildet; die Spalten von \(A_i\) sind die farbigen Kantenvektoren.

3D: Tetraeder

Das Referenztetraeder \(T\) hat die Ecken \((0,0,0)\), \((1,0,0)\), \((0,1,0)\) und \((0,0,1)\). Ein Punkt in \(T\) wird mit \(\vec{t}=(\xi,\eta,\zeta)^T\) bezeichnet. Mit den Kantenvektoren

\[ \vec{a} = \vec{p}_{i,2}-\vec{p}_{i,1},\qquad \vec{b} = \vec{p}_{i,3}-\vec{p}_{i,1},\qquad \vec{c} = \vec{p}_{i,4}-\vec{p}_{i,1} \]

gilt Folgendes.

Abbildung.

\[\begin{split} \vec{x} = \sigma_i(\vec{t}\,) = \sigma_i(\xi,\eta,\zeta) = \vec{p}_{i,1} + A_i\cdot\begin{pmatrix}\xi\\\eta\\\zeta\end{pmatrix}, \qquad A_i = (\vec{a},\ \vec{b},\ \vec{c}), \end{split}\]

bzw. \(\sigma_i(\xi,\eta,\zeta) = (1-\xi-\eta-\zeta)\,\vec{p}_{i,1} + \xi\,\vec{p}_{i,2} + \eta\,\vec{p}_{i,3} + \zeta\,\vec{p}_{i,4}\).

Jacobi-Determinante (Spatprodukt).

\[ J_i = \det A_i = \vec{a}\cdot(\vec{b}\times\vec{c}) = \pm 6\,|T_i| , \qquad dx\,dy\,dz = |J_i|\,d\xi\,d\eta\,d\zeta . \]

\(J_i\) ist genau dann positiv, wenn \((\vec{a},\vec{b},\vec{c})\) ein Rechtssystem bilden.

Inverse Jacobi-Matrix. Die Zeilen von \(A_i^{-1}\) sind (bis auf den Faktor \(1/J_i\)) Kreuzprodukte der Kantenvektoren:

\[\begin{split} A_i^{-1} = \frac{1}{J_i} \begin{pmatrix} (\vec{b}\times\vec{c})^T\\ (\vec{c}\times\vec{a})^T\\ (\vec{a}\times\vec{b})^T \end{pmatrix}. \end{split}\]

Probe: \((\vec{b}\times\vec{c})\cdot\vec{a} = J_i\) und \((\vec{b}\times\vec{c})\cdot\vec{b} = (\vec{b}\times\vec{c})\cdot\vec{c}=0\). Geometrisch sind die Zeilen also Normalenvektoren der Seitenflächen durch \(\vec{p}_{i,1}\).

Rücktransformation.

\[\begin{split} \begin{pmatrix}\xi\\\eta\\\zeta\end{pmatrix} = \sigma_i^{-1}(\vec{x}) = A_i^{-1}\cdot\bigl(\vec{x}-\vec{p}_{i,1}\bigr), \qquad \begin{aligned} \xi &= \frac{(\vec{b}\times\vec{c})\cdot(\vec{x}-\vec{p}_{i,1})}{J_i},\\ \eta &= \frac{(\vec{c}\times\vec{a})\cdot(\vec{x}-\vec{p}_{i,1})}{J_i},\\ \zeta &= \frac{(\vec{a}\times\vec{b})\cdot(\vec{x}-\vec{p}_{i,1})}{J_i}. \end{aligned} \end{split}\]

Nach Baryzentrische Koordinaten ist \(\xi=\lambda_2\), \(\eta=\lambda_3\) und \(\zeta=\lambda_4\). Zum Beispiel ist \((\vec{b}\times\vec{c})\cdot(\vec{x}-\vec{p}_{i,1}) = \det(\vec{x}-\vec{p}_{i,1},\ \vec{b},\ \vec{c})\) das Sechsfache des orientierten Volumens des Teiltetraeders \((\vec{p}_{i,1},\vec{x},\vec{p}_{i,3},\vec{p}_{i,4})\), in dem \(\vec{p}_{i,2}\) durch \(\vec{x}\) ersetzt ist. Geteilt durch \(J_i = 6\,|T_i|\) ergibt sich der Volumenanteil.

Beispiel. \(\vec{p}_{i,1}=(1,1,0)^T\), \(\vec{p}_{i,2}=(3,1,0)^T\), \(\vec{p}_{i,3}=(1.5,3,0)^T\), \(\vec{p}_{i,4}=(1.5,1.5,2)^T\):

\[\begin{split} A_i = \begin{pmatrix} 2 & 0.5 & 0.5\\ 0 & 2 & 0.5\\ 0 & 0 & 2\end{pmatrix},\qquad J_i = 8,\qquad |T_i| = \tfrac{8}{6} = \tfrac43, \end{split}\]

\[\begin{split} A_i^{-1} = \begin{pmatrix} \tfrac12 & -\tfrac18 & -\tfrac{3}{32}\\ 0 & \tfrac12 & -\tfrac18\\ 0 & 0 & \tfrac12\end{pmatrix} = \begin{pmatrix} 0.5 & -0.125 & -0.09375\\ 0 & 0.5 & -0.125\\ 0 & 0 & 0.5\end{pmatrix}. \end{split}\]

Probe: \(\sigma_i^{-1}(\vec{p}_{i,4}) = A_i^{-1}\cdot(0.5,\,0.5,\,2)^T = (0,0,1)^T\). Der Schwerpunkt \(\tfrac14\sum_k\vec{p}_{i,k} = (1.75,\,1.625,\,0.5)^T\) wird auf \((\tfrac14,\tfrac14,\tfrac14)^T\) abgebildet.

Abb. 9.11 Affine Abbildung zwischen einem beliebigen Tetraeder \(T_i\) und dem Einheitstetraeder \(T\) (verdeckte Kanten gestrichelt). Der Schwerpunkt wird auf den Schwerpunkt abgebildet.

Numerische Kontrolle

Die Beispiele lassen sich für alle Dimensionen mit derselben Funktion nachrechnen:

import numpy as np
from math import factorial

def affine_map(pts):
    """Jacobi-Matrix A, sigma und sigma^{-1} für ein Simplex mit Ecken pts (Zeilen)."""
    pts = np.asarray(pts, dtype=float)
    p1 = pts[0]
    A = (pts[1:] - p1).T                 # Spalten = Kantenvektoren
    Ainv = np.linalg.inv(A)
    sigma = lambda t: p1 + A @ t         # Referenz -> Element
    sigma_inv = lambda x: Ainv @ (x - p1)  # Element -> Referenz
    return A, Ainv, sigma, sigma_inv

for pts in ([[1], [4]],
            [[1, 1], [4, 2], [2, 4]],
            [[1, 1, 0], [3, 1, 0], [1.5, 3, 0], [1.5, 1.5, 2]]):
    A, Ainv, sigma, sigma_inv = affine_map(pts)
    d = len(pts) - 1
    J = np.linalg.det(A)
    centroid = np.mean(pts, axis=0)
    print(f"d={d}: J = {J:.4g}, |T_i| = {abs(J)/factorial(d):.4g}, "
          f"sigma^-1(Schwerpunkt) = {sigma_inv(centroid)}")

Zusammenfassung

	1D: Intervall	2D: Dreieck	3D: Tetraeder
\(A_i\)	\(b-a\)	\((\vec{p}_{i,2}-\vec{p}_{i,1},\ \vec{p}_{i,3}-\vec{p}_{i,1})\)	\((\vec{a},\ \vec{b},\ \vec{c})\)
\(J_i=\det A_i\)	\(b-a\)	\(\pm 2\vert T_i\vert\)	\(\vec{a}\cdot(\vec{b}\times\vec{c})=\pm 6\vert T_i\vert\)
\(\sigma_i(\vec{t}\,)\)	\(a+(b-a)\xi\)	\(\vec{p}_{i,1}+A_i\cdot\vec{t}\,\)	\(\vec{p}_{i,1}+A_i\cdot\vec{t}\,\)
\(\sigma_i^{-1}(\vec{x})\)	\(\frac{x-a}{b-a}\)	\(A_i^{-1}\cdot(\vec{x}-\vec{p}_{i,1})\)	\(A_i^{-1}\cdot(\vec{x}-\vec{p}_{i,1})\)
\(\vert T\vert\)	\(1\)	\(\tfrac12\)	\(\tfrac16\)

Für alle Fälle gilt

\[\begin{split} \begin{gathered} \int_{T_i} f(\vec{x})\,d\vec{x} = \int_{T} f\bigl(\sigma_i(\vec{t}\,)\bigr)\,|J_i|\,d\vec{t},\qquad \nabla\varphi(\vec{x}) = {A_i^{-1}}^T\cdot\nabla\tilde{\varphi}\bigl(\sigma_i^{-1}(\vec{x})\bigr),\\ \int_{T_i} \nabla\varphi_j(\vec{x})\cdot\nabla\varphi_k(\vec{x})\,d\vec{x} = \int_{T} \nabla\tilde{\varphi}_j(\vec{t}\,)^T\cdot C_i\cdot\nabla\tilde{\varphi}_k(\vec{t}\,)\;|J_i|\,d\vec{t}, \qquad C_i = A_i^{-1}\cdot{A_i^{-1}}^T . \end{gathered} \end{split}\]

Da \(\sigma_i\) affin ist, sind \(A_i\), \(J_i\) und \(C_i\) auf ganz \(T_i\) konstant. Deshalb genügt es, Basisfunktionen und Quadraturregeln nur einmal auf \(T\) zu definieren.