Praktikum Wissenschaftliches Rechnen: Elektrische Impedanztomographie

Wintersemester 1998/99

Dr. A. Klawonn^{1}  Dr. F. Wübbeling^{2}

12. November 1998

1e-mail: klawonn@math.uni-muenster.de 2e-mail: wuebbel@math.uni-muenster.de

2

Inhaltsverzeichnis

1 Finite Elemente 1 1.1 Variationsformulierungen elliptischer Differentialgleichungen . . . . . . 1 1.2 Das Lemma von Lax-Milgram . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Das Galerkin-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Das einzelne finite Element . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.5 Forderungen an die Triangulierung . . . . . . . . . . . . . . . . . . . . 6 1.6 Beispiele Finiter Elemente . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.6.1 Lineare Dreiecks-Elemente (Lagrange ) . . . . . . . . . . . . . . 7 1.7 Zusammensetzen Finiter Elemente . . . . . . . . . . . . . . . . . . . . . 7 1.8 Implementierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.8.1 Gittererzeugung . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.8.2 Aufstellen des Gleichungssystems . . . . . . . . . . . . . . . . . 10

Bibliography 19

ii

Kapitel 1

Finite Elemente

In diesem Kapitel soll ein kurzer Überblick zur Methode der finiten Elemente gegeben werden. Dabei liegt der Schwerpunkt auf der Implementierung. Für eine detaillierte theoretische Behandlung verweisen wir auf die Lehrbücher von Braess [2] und Brenner/Scott [3].

1.1 Variationsformulierungen elliptischer Differentialgleichungen

Wir beschränken uns in dieser Einführung auf gleichmäßig elliptische partielle Differentialgleichungen zweiter Ordnung. Weiterhin nehmen wir an, daß der Leser mit den Grundlagen der Sobolevräume vertraut ist, wie sie zum Beispiel in [5] vermittelt werden, s.a. [2, 3].

Sei Ω ⊂ R^{d}, d = 2, 3 ein Gebiet, dann betrachten wir die partielle Differentialgleichung zweiter Ordnung

−div(a∇u) + b u = f in Ω
u = g_{0} auf Γ_{0}
\frac{∂u}{∂ν} + c u = g_{1} auf Γ_{1},

(1.1)

wobei Γ_{0}, Γ_{1} ⊂ ∂Ω mit Γ_{0} ∪ Γ_{1} = ∂Ω und Γ_{0} ∩ Γ_{1} = ∅ gilt. Desweiteren sind a = a(x), b = b(x) und c = c(x) gegebene Funktionen der Ortsvariablen x. Gesucht ist die unbekannte Funktion u = u(x).

Bei den Randbedingungen können auch die Extremfälle Γ_{0} = ∂Ω und Γ_{1} = ∂Ω auftreten.

Wir benutzen den Sobolevraum H^{1}(Ω) = \overline{C^{∞}(Ω)}^{∥·∥_{H^{1}(Ω)}}, wobei der Abschluß bezüglich der Norm ∥u∥_{H^{1}(Ω)} := \sqrt{|u|_{H^{1}(Ω)}^{2} + ∥u∥_{L_{2}(Ω)}^{2}} mit

|u|_{H^{1}(Ω)}^{2} := ∫_{Ω}(∇u)^{2}dx und ∥u∥_{L_{2}(Ω)}^{2} := ∫_{Ω} u^{2}dx

gebildet wird. Desweiteren führen wir den Soblevraum H_{0}^{1}(Ω) := \overline{C_{0}^{∞}(Ω)}^{∥·∥_{H^{1}(Ω)}} ein.

Da wir uns mit L_{2}−Funktionen beschäftigen, müssen wir noch erklären, was wir unter Randwerten solcher Funktionen verstehen. Es gilt folgender Satz

Satz 1.1 (Spursatz) Sei Ω ein Lipschitz–Gebiet. Dann existiert eine stetige, lineare Abbildung

γ : H^{1}(Ω) → L_{2}(∂Ω)

mit

γ u(x) = u(x) ∀x ∈ ∂Ω,

für alle u ∈ C^{1}(Ω̄).

Bemerkung 1.1 Anschaulich bedeutet der Spursatz, dass γu die Einschränkung von u auf den Rand ∂Ω ist. Die Auswertung einer H^{1}−Funktion an einem Punkt ist nicht immer sinnvoll, aus dem Spursatz folgt aber immerhin, dass die Restriktion einer H^{1}−Funktion immerhin noch eine L_{2}−Funktion ist.

Beschränken wir uns auf Lipschitz-Gebiete, so gilt

H_{0}^{1}(Ω) = {v ∈ H^{1}(Ω) : v_{|∂Ω} = 0},

d.h. H_{0}^{1}(Ω) besteht aus genau den H^{1}−Funktionen, die auf dem Rand von Ω den Wert Null haben. Sei Ω ein Lipschitz-Gebiet, dann definieren wir schließlich noch die Räume

H_{Γ_{0},g_{0}}^{1}(Ω) := {v ∈ H^{1}(Ω) : v_{|Γ_{0}} = g_{0}},
H_{Γ_{0},0}^{1}(Ω) := {v ∈ H^{1}(Ω) : v_{|Γ_{0}} = 0}.

Damit die Definition von H_{Γ_{0},g_{0}}^{1}(Ω) wohldefiniert ist, nehmen wir der Einfachheit halber g_{0} ∈ H^{1}(Γ_{0}) an. Für genauere Aussagen benötigen wir die Theorie der Sobolevräume H^{k} mit reellem Exponenten k, s. zum Beispiel Lions/Magenes [8], Nečas [9], Wloka [10], Hackbusch [7]. Es gilt folgender

Satz 1.2 Die Sobolevräume H^{1}(Ω), H_{0}^{1}(Ω) und H_{Γ_{0},0}^{1}(Ω) sind jeweils zusammen mit dem inneren Produkt (u, v)_{H^{1}(Ω)} := (∇u, ∇v)_{L_{2}(Ω)} + (u, v)_{L_{2}(Ω)}, (u, v)_{L_{2}(Ω)} := ∫_{Ω} uvdx Hilberträume.

Bemerkung 1.2 Man beachte, dass H_{Γ_{0},g_{0}}^{1}(Ω) noch nicht einmal ein linearer Vektorraum ist. Begründung!

Unter Benutzung des Gauß’schen Integralsatzes und des üblichen Dichteschluss können wir eine Variationsformulierung des Randwertproblems (1.1) herleiten. Das Problem lautet dann:

Finde ein u ∈ H_{Γ_{0},g_{0}}^{1}(Ω), so dass

A(u, v) = F(v) ∀v ∈ H_{Γ_{0},0}^{1}(Ω),  (1.2)

wobei

A(u, v) = ∫_{Ω} a∇u∇v + buvdx + ∫_{Γ_{1}} a c u vds,
=: A_{Ω}(u, v) + A_{Γ_{1}}(u, v),
F(v) = ∫_{Ω} f vdx + ∫_{Γ_{1}} a g_{1} vds,
=: F_{Ω}(v) + F_{Γ_{1}}(v).

Bei dieser Formulierung müssen wir noch mit zwei verschiedenen Räumen arbeiten, nämlich H_{Γ_{0},0}^{1}(Ω) und H_{Γ_{0},g_{0}}^{1}(Ω). Kennen wir jedoch eine Funktion u_{g} ∈ H_{Γ_{0},g_{0}}^{1}(Ω), so können wir H_{Γ_{0},g_{0}}^{1}(Ω) wie folgt schreiben

H_{Γ_{0},g_{0}}^{1}(Ω) = {u ∈ H^{1}(Ω) : u = u_{g} + v, v ∈ H_{Γ_{0},g_{0}}^{1}(Ω)}.

Damit ergibt sich folgende zu (1.2) äquivalente Variationsformulierung:

Finde u_{0} ∈ H_{Γ_{0},0}^{1}(Ω), so dass

A(u_{0} + u_{g}, v) = F(v) ∀v ∈ H_{Γ_{0},0}^{1}(Ω).  (1.3)

Die Lösung von (1.2) ergibt sich dann als u = u_{0} + u_{g}. Die Formulierung (1.3) führt also auf das Ersatzproblem

Finde u_{0} ∈ H_{Γ_{0},0}^{1}(Ω), so dass

A(u_{0}, v) = F̂(v) ∀v ∈ H_{Γ_{0},0}^{1}(Ω),  (1.4)

mit F̂(v) := F(v) − A(u_{g}, v).

Um das ursprüngliche Randwertproblem (1.1) zu lösen, muß also zunächst eine beliebige Funktion u_{g} ∈ H_{Γ_{0},g_{0}}^{1}(Ω) konstruiert werden. Die Lösung ergibt sich dann mit Hilfe von (1.4).

1.2 Das Lemma von Lax-Milgram

In dieser Sektion betrachten wir die Eindeutigkeit und Lösbarkeit von Variationsgleichungen in Hilberträumen. Wir werden sehen, dass unter gewissen, auch in der Praxis sinnvollen, Voraussetzungen die Variationsaufgabe (1.4) ein gut gestelltes Problem im Sinne von Hadamard ist.

Sei V ein Hilbertraum mit innerem Produkt (·, ·) und der dadurch induzierten, zugehörigen Norm ∥u∥ := \sqrt{(u, u)}, weiterhin seien A(·, ·) : V × V → R eine stetige Bilinearform auf V und F(·) : V → R ein stetiges lineares Funktional auf V. Wir betrachten folgendes lineares Variationsproblem:

Finde ein u ∈ V, so dass

A(u, v) = F(v) ∀v ∈ V.  (1.5)

Das folgende Lemma von P.D. Lax und A. Milgram ist eine einfache, aber wirkungsvolle Verallgemeinerung des Riesz’schen Darstellungssatzes, für einen Beweis vgl. [2, 3].

Lemma 1.1 (Lax-Milgram) Die Bilinearform A(·, ·) sei stetig, d.h.

∃M ≡ const > 0 |A(u, v)| ≤ M ∥u∥∥v∥ ∀u, v ∈ V,

und V −elliptisch, d.h.

∃α ≡ const > 0 A(u, u) ≥ α ∥u∥^{2} ∀u ∈ V.

Dann hat die Variationsgleichung (1.5) eine eindeutige Lösung u und es gilt

∥u∥ ≤ \frac{1}{α}∥F∥.

Aus der letzten Ungleichung im Lemma von Lax-Milgram folgt die stetige Abhängigkeit der Lösung von den rechten Seiten. Somit liegt ein gut gestelltes Problem im Sinne von Hadamard vor.

Nach Satz 1.2 sind die Räume H^{1}(Ω), H_{0}^{1}(Ω) und H_{Γ_{0},0}^{1}(Ω) Hilberträume, jedoch nicht der Raum H_{Γ_{0},g_{0}}^{1}(Ω) !

Aufgabe 1.1 Man wende das Lemma von Lax–Milgram auf die Variationsprobleme der Sektion 1.1 an. Dabei soll der Leser sich überlegen, welche Normen und inneren Produkte benutzt werden müssen.

1.3 Das Galerkin-Verfahren

In diesem Abschnitt betrachten wir ein Näherungsverfahren zur Lösung von Variationsgleichungen: das Galerkin–Verfahren. Im Gegensatz zum Verfahren der finiten Differenzen, bei dem der Differentialoperator diskretisiert wird, betrachten wir beim Galerkin–Verfahren den kontinuierlichen Operator bzw. die zugehörige Bilinearform und diskretisieren den Testraum V.

Wir stellen im Folgenden dieselben Anforderungen an die Bilinearform A(·, ·), die Linearform F(·) und den Raum V. Mit u ∈ V bezeichnen wir die Lösung der Variationsgleichung (1.5).

Sei V_{h} ein endlichdimensionaler Unterraum von V. Dann heißt u_{h} Galerkin–Näherung von u, falls

A(u_{h}, v) = F(v) ∀v ∈ V_{h}.  (1.6)

Ist (v_{i})_{i=1,...,M} eine Basis des Teilraums V_{h}, so existiert eine Entwicklung von u_{h} in dieser Basis, d.h. es existieren reelle Koeffizienten (d_{k})_{k=1,...,M} mit

u_{h} = ∑_{i=1}^{M} d_{i}v_{i}.

Statt das Problem (1.6) zu betrachten, genügt es

A(u_{h}, v_{i}) = F(v_{i}) ∀v_{i}, i = 1, …, M,  (1.7)

zu lösen.

Setzen wir die Basisentwicklung u_{h} = ∑_{i=1}^{M} d_{i}v_{i} in die Variationsformulierung (1.7) ein, so erhalten wir das lineare Gleichungssystem

K d = F  (1.8)

mit K = (A(v_{i}, v_{j}))_{i,j=1,...,M}, d = (d_{i})_{i=1,...,M} und F = (F (v_{i}))_{i=1,...,M}.

Beim Galerkin–Verfahren löst man die Variationsgleichung (1.5) durch Einschränkung auf einen endlichdimensionalen Teilraum und das Aufstellen und Lösen eines linearen Gleichungssystems.

Der nächste Satz zeigt, dass die Galerkin–Näherung die Bestapproximation aus V_{h} an u ∈ V ist.

Satz 1.3 Die Voraussetzungen von Lemma 1.1 seien erfüllt, dann existiert u_{h} und ist eindeutig bestimmt. Weiterhin gelten die Abschätzungen

∥u_{h}∥ ≤ \frac{1}{α} ∥F∥, ∥u − u_{h}∥ ≤ \frac{M}{α} min_{v∈V_{h}} ∥u − v∥  (1.9)

Das Verfahren der finiten Elemente ist ein Galerkin–Verfahren, bei dem die Basisfunktionen v_{i}, i = 1, …, N mit Hilfe einer Zerlegung des Gebiets Ω in kleinere Teilgebiete mit dem Durchmesser h konstruiert werden. In diesen kleinen Teilgebieten werden die Funktionen aus V durch einfache Funktionen mit kleinem Träger, in der Regel Polynome niedrigen Grades, approximiert.

1.4 Das einzelne finite Element

Wir benutzen hier die Definition eines finiten Elements nach Ciarlet [4].

Definition 1.1 Ein Finites Element ist ein Tripel (E, P, N), wobei (i) E ⊂ R^{d} ein Lipschitz-Gebiet,

(ii) P ein endlichdimensionaler Teilraum von C(\overline{E}), dem Raum der stetigen Funktionen auf \overline{E} und

(iii) N = {N_{1}, N_{2}, …, N_{m}} eine Basis für P′, den Dualraum von P,

ist.

Bemerkung 1.3 (i) E ist in der Regel ein einfaches Gebiet, z.B. ein Dreieck (bzw. Tetraeder) oder ein Viereck (bzw. Quader). Häufig bezeichnet man auch E selbst als Finites Element oder einfach als Element.

(ii) P besteht in der Regel aus Polynomen niedrigen Grades. Die Dimension m von P ist üblicherweise klein, etwa 3 − 20.

(iii) N besteht meist aus Funktionsauswertungen bzw. Auswertung von Ableitungen in bestimmten Punkten von \overline{E}, den Knoten. Die Funktionale in N nennt man auch Freiheitsgrade des Elements.

Definition 1.2 Sei (E, P, N) ein Finites Element und sei (ϕ_{i})_{i=1,...,m} eine Basis von P, die dual ist zu N, d.h. N_{i}(ϕ_{j}) = δ_{ij}, i, j = 1, …, m. Die Basis (ϕ_{i})_{i=1,...,m} heißt Knotenbasis oder nodale Basis.

1.5 Forderungen an die Triangulierung

Wir nehmen im folgenden der Einfachheit halber an, daß unser Gebiet Ω polygonal ist.

Definition 1.3 Eine Zerlegung oder Triangulierung τ = {E_{1}, …, E_{L}} von Ω heißt zulässig, wenn folgendes gilt:

(i) Ω̄ = ⋃_{l=1}^{L} E_{l},

(ii) Besteht E_{l} ∩ E_{k} aus genau einem Punkt, so ist dieser Eckpunkt sowohl von E_{l} als auch von E_{k},

(iii) Besteht E_{l}∩E_{k}, k ≠ l aus mehr als einem Punkt, so handelt es sich um eine Kante oder Seitenfläche sowohl von E_{l} als auch von E_{k}.

1.6 Beispiele Finiter Elemente

Zunächst beschränken wir uns auf zweidimensionale Elemente. Zur Beschreibung der verschiedenen Elemente geben wir jeweils das zugehörige Tripel (E, P, N) an.

1.6.1 Lineare Dreiecks-Elemente (Lagrange )

Das Element-Gebiet ist hier ein Dreieck mit den Ecken a_{i} = (x_{i}, y_{i}), i = 1, 2, 3. Wir wählen P = P_{1}, d.h. den Raum der linearen Polynome, also m = dimP = 3 und wir setzen N = {N_{1}, N_{2}, N_{3}} mit N_{i}(p) = p(a_{i}), i = 1, 2, 3.

Für das Einheitsdreieck ergibt sich als Knotenbasis (ϕ_{i})_{i=1,2,3}

ϕ_{1}(x, y) := 1 − x − y, ϕ_{2}(x, y) := x, ϕ_{3}(x, y) := y.

1.7 Zusammensetzen Finiter Elemente

Sei Ω ⊂ R^{d} ein Lipschitz-Gebiet und seien (E_{l}, P_{l}, N_{l}), l = 1, …, L Finite Elemente mit

\overline{Ω} = ⋃_{l=1}^{L} \overline{E}_{l}

und E_{l} ∩ E_{k} = ∅ für l ≠ k.

Definition 1.4 Eine Familie Finiter Elemente (E_{l}, P_{l}, N_{l}), l = 1, …, L heißt affin, wenn jedes Element der Familie affines Bild eines einzigen Elementes ist. Dieses Element wird auch als Referenzelement bezeichnet.

Beispiel 1.1 Bei linearen Dreieckselementen wählt man üblicherweise das Einheitsdreieck als Referenzelement. Die affine lineare Abbildung lautet dann

T : Ê → E
T x̂ := Bx̂ + b

mit

B = [ \begin{matrix}x_{2} − x_{1} & x_{3} − x_{1} \\ y_{2} − y_{1} & y_{3} − y_{1} \end{matrix} ] , b = [ \begin{matrix}x_{1} \\ y_{1} \end{matrix} ] ,

wobei Ê das Einheisdreieck bezeichnet und das jeweilige Element E durch die Ecken a_{i} = (x_{i}, y_{i}), i = 1, 2, 3 gegeben ist.

Definition 1.5

(i) Eine Familie (E_{i}, P_{i}, N_{i})_{i=1,...,L} heisst verträglich, wenn für jeden Knoten a, der sowohl zu E_{i} als auch zu E_{k} gehört gilt: span{N_{ai}} = span{N_{ak}}, wobei N_{ai} := {N_{j}^{i} ∈ N_{i} : N_{j}^{i} gehört zu Knoten a}.

(ii) Ist eine Familie verträglich, so lassen sich zu jedem Knoten a dieser Familie Freiheitsgrade φ_{1}, …, φ_{t} finden, so dass span{φ_{1}, …, φ_{t}} = span{N_{ai}} für jedes Element E_{i}, zu dem der Knoten a gehört. Die Menge der (φ_{j})_{j=1,...,M}, bildet die Menge der globalen Freiheitsgrade.

1.8 Implementierung

In diesem Abschnitt soll die übliche Vorgehensweise bei der Programmierung des Verfahrens der Finiten Elemente beschrieben werden. Wir verwenden dabei die Bezeichungen der vorhergehenden Sektionen. Das Lösen einer Differentialgleichung mit Hilfe der Finiten Elemente kann in drei Abschnitte gegliedert werden:

(i) Gittererzeugung

(ii) Aufstellen des Gleichungssystems (Steifigkeitsmatrix, Lastvektor)

(iii) Lösen des Gleichungssystems und evtl. Bearbeitung der Lösung (Visualisierung, etc.)

1.8.1 Gittererzeugung

Es gibt viele Artikel, Bücher und Programme zur Gittererzeugung für geometrisch komplexe Gebiete, s. z.B. [6]. Wir beschränken uns der Einfachheit halber jedoch auf leicht zu zerlegende Gebiete, z.B. Kreis oder Rechteck, bzw. Gebiete, die sich daraus zusammensetzen.

Unabhängig von der Geometrie des Gebiets benötigen wir gewisse Daten der Triangulierung. Wir folgen im Wesentlichen den Bezeichnungen in Axelsson/Barker [1][Kap. 5]. Sei Ω ein Gebiet und es sei ein Gitter mit L Elementen und M Knoten gegeben. Weiterhin sei {Φ_{1}, …, Φ_{M}} die zugehörige Knotenbasis, vgl. Definition 1.2. Desweiteren unterscheiden wir zwischen Knoten, die im Innern des Gebiets oder auf einem Teil des Randes mit vorgeschriebenen Neumann-Randwerten liegen und Knoten, die auf einem Randstück mit vorgeschriebenen Dirichlet-Randdaten liegen. Dies ist sinnvoll, da die Lösung in den Knoten mit Dirichletbedingungen schon bekannt sind, also anders behandelt werden müssen als ”echte“ Unbekannte. Im Folgenden bezeichnen wir die Anzahl der Knoten auf Γ_{0} mit M_{0}, die Indexmenge der Knotennummern von Knoten auf Γ_{0} mit J_{0}. Die Anzahl der Knoten in Ω ∪ Γ_{1} ist N und die Indexmenge aller Knotennummern von Knoten in Ω ∪ Γ_{1} wird mit J bezeichnet.

Zum Aufstellen des Gleichungssystems müssen alle Knoten durchlaufen werden. Dabei hat es sich als sinnvoll herausgestellt, elementweise vorzugehen, d.h. es werden nacheinander alle Elemente durchlaufen und pro Element werden die zugehörigen Knoten abgearbeitet. Natürlich werden bei dieser Vorgehensweise die meisten Knoten mehrmals besucht. Um dieses Verfahren zu implementieren, benötigen wir also eine globale Knoten- und eine globale Elementnummerierung, sowie eine lokale Knotennummerierung pro Element. Die global Knotennummerierung kann beliebig erfolgen, sie muss nur die Bedingung erfüllen, dass jedem Knoten genau eine (globale) Nummer zugeordnet wird. Analog müssen die Elemente durchnummeriert werden. Innerhalb jedes Elements führen wir noch eine lokale Knotennummerierung ein, die - für jedes Element gleich - die Knoten, welche zum Element gehören, (lokal) nummeriert. Bei affinen Finite-Element-Familien lässt sich die lokale Knotennummerierung bequem über das Referenzelement definieren. Zur Berechnung der Steifigkeitsmatrix und des Lastvektors brauchen wir auch die Koordinaten der jeweiligen Knoten, schließlich müssen auch noch die Daten der Randknoten erfasst werden. Diese Daten werden alle in verschiedenen Listen gespeichert, die wir hier der Einfachheit halber als Arrays darstellen:

ELNODE L × m Array, das die globalen Knotennummern

für alle Elemente enthält

NODECO M × 2 Array, das die Koordinaten aller Knoten enthält
BONODE M_{0} × 1 Array, das die Indexmenge J_{0} enthält
EDGE  L_{1} × 2 enthält die Daten der Elemente, die auf Γ_{1} liegen.

Die Arrays sind dabei wie folgt aufgebaut:

ELNODE  NODECO
Element-Nr. Element-Knoten  Knoten-Nr. Knoten-Koordinaten
1  E_{1}^{1}, …, E_{m}^{1}  1  (x_{1}, y_{1})
⋮  ⋮  ⋮  ⋮
L  E_{1}^{L}, …, E_{m}^{L}  M  (x_{M}, y_{M})

Dabei enthält E_{i}^{l}, i = 1, …, m, l = 1, …, L die globale Knotennummer des i−ten lokalen Knotens im l−ten Element.

Die Randdaten werden in den folgenden beiden Listen gespeichert:

BONODE  EDGE
Randknoten-Nr. Knoten-Nr.(global)  Randelement-Nr. Knoten (lokal)
1  i_{1}  1  RE_{1}^{First}
⋮  ⋮  ⋮  ⋮
M_{0}  i_{M_{0}}  L_{1}  RE_{L_{1}}^{First}

Die Reihenfolge der Randknoten in BONODE und die der Randelemente in ED-GE ist willkürlich. Die Liste BONODE enthält zu jeder Randknotennummer die entsprechende globale Knotennummer. In EDGE findet sich zu jeder Randelementnummer die zugehörige lokale Knotennummer des ersten Randknotens des jeweiligen Elements. Mit ”ersten“ ist gemeint, dass alle weiteren Randknoten dieses Elements im mathematisch positiven Drehsinn nach diesem Knoten kommen. Wir weisen darauf hin, dass eine Ecke des Gebiets Ω zu Elementen führen kann, die mehr als einen Elementrand auf Γ_{1} haben.

Beispiel 1.2 Vorgehensweise zur Triangulierung eines Kreises wie auf Aufgabenblatt 1 beschrieben

1.8.2 Aufstellen des Gleichungssystems

Wir betrachten nun das Variationsproblem (1.6)

A(u_{h}, v_{h}) = F(v_{h}) ∀v_{h} ∈ V_{h}.  (1.10)

Sei eine affine verträgliche Familie Finiter Elemente (E_{i}, P_{i}, N_{i})_{i=1,...,L} gegeben. Die zugehörigen globalen Freiheitsgrade bezeichnen wir mit (φ_{j})_{j=1,...,M}, die globalen Knoten– Basisfunktionen mit (ϕ_{j})_{j=1,...,M}. Somit gilt V_{h} = span_{j=1}^{M}{ϕ_{j}}.

Um unser Programm möglichst modular zu halten, werden wir die Implementierung der Randbedingungen getrennt behandeln und uns zunächst nur mit den Teilen des Programms beschäftigen, bei denen über Ω integriert wird.

Sei (E, P, N) ein beliebiges Element der o.a. Familie und (Φ_{E,j})_{j=1,...,m} die lokale Knotenbasis. Dann gilt für u_{h}, v_{h} ∈ V_{h} eingeschränkt auf das Element E

u_{h} = ∑_{j=1}^{m} u_{E,j}Φ_{E,j}, v_{h} = ∑_{j=1}^{m} v_{E,j}Φ_{E,j},

wobei mit u_{E,j}, v_{E,j} die Werte von u_{h}, v_{h} in den Knoten des Elements E bezeichnet werden.

Wir betrachten zunächst nur Anteile der Variationsgleichung ohne Randbedingungen, d.h. A_{Ω}(·, ·) und F_{Ω}(·). Beschränken wir weiterhin den Integrationsbereich auf E und benutzen die Entwicklung von u_{h}, v_{h} in die lokale Knotenbasis, so erhalten wir

A_{E}(u_{h}, v_{h}) = ∑_{i,j=1}^{m} u_{E,i}v_{E,j}A_{E}(Φ_{E,i}, Φ_{E,j})
= v_{E,h}^{t}K_{E}u_{E,h},

wobei A_{E}(·, ·) aus A_{Ω}(·, ·) durch Einschränkung des Integrationsgebiets auf E entsteht. Desweiteren benutzen wir

u_{E,h} := [ u_{E,1} ; ⋮ ; u_{E,m} ] , v_{E,h} := [ v_{E,1} ; ⋮ ; v_{E,m} ] , K_{E} := [ \begin{matrix}A_{E}(Φ_{E,1}, Φ_{E,1}) & … & A_{E}(Φ_{E,1}, Φ_{E,m}) \\ ⋮ & ... & ⋮ \\ A_{E}(Φ_{E,m}, Φ_{E,1}) & … & A_{E}(Φ_{E,m}, Φ_{E,m} \end{matrix} ] .

Analog ergibt sich

F_{E}(v_{h}) = v_{E,h}^{t}F_{E}u_{E,h},

mit

F_{E} := [ F_{E}(Φ_{E,1}) ; ⋮ ; F_{E}(Φ_{E,m} ] .

Die Matrix K_{E} und der Vektor F_{E} heißen Element–Matrix und –Vektor. Da die betrachtete Familie verträglich ist, gilt

A_{Ω}(u_{h}, v_{h}) = ∑_{l=1}^{L} v_{E_{l},h}^{t}K_{E_{l}}u_{E_{l},h}, F_{Ω} = ∑_{l=1}^{L} v_{E_{l},h}^{t}F_{E_{l}}.

Für jedes Element E_{l}, l = 1, …, L lässt sich eine Indexmenge e_{i}^{l}, i = 1, …, m finden, so dass der i−te (lokale) Knoten des Elements E_{l} der e_{i}^{l}−te Knoten in der globalen Knotennummerierung ist. Diese Indexmenge ist in der Liste ELNODE des Abschnitts 1.8.1 realisiert. Wir führen nun für jedes Element der Triangulierung Restriktionsmatrizen R_{l} ein. Dabei handelt es sich um Boole’sche Matrizen, die wie folgt definiert sind:

R_{l} := (R_{ij}^{l})_{{i=1,...,m \atop j=1,...,M}} , R_{ij}^{l} := δ_{e_{i}^{l}j}.

Verwenden wir die Bezeichnungen aus dem Abschnitt (1.8.1), so können wir die Werte der Funktionen u_{h}, v_{h} in den Knoten zu Vektoren u, v zusammenfassen:

u := [ u_{1} ; ⋮ ; u_{M} ] , v := [ v_{1} ; ⋮ ; v_{M} ] ,

Dabei sind die Einträge u_{i}, v_{i} der Vektoren u, v wie folgt definiert

u_{e_{i}^{l}} := u_{E_{l},i}, v_{e_{i}^{l}} := v_{E_{l},i}.

Offenbar gilt nun

A_{Ω}(u_{h}, v_{h}) = v^{t}Ku, F_{Ω}(v_{h}) = v^{t}F,

wobei

K := ∑_{l=1}^{L} R_{l}^{t}K_{E_{l}}R_{l}, F := ∑_{l=1}^{L} R_{l}^{t}F_{E_{l}}.

Die Matrix K heißt Steifigkeits– oder Systemmatrix, der Vektor F Last– oder Systemvektor. Man beachte, dass wir die Randbedingungen noch nicht berücksichtigt haben. Wir werden jedoch die Namen Steifigkeits-/Systemmatrix bzw. Last-/Systemvektor für beide Fälle benutzen.

Bei der Einarbeitung der Randbedingungen kann ähnlich wie in Abschnitt 1.1 vorgegangen werden, vgl. (1.4).

Das Hinzufügen der Neumann-Randbedingungen ist prinzipiell einfach und lässt sich durch einfaches Addieren der Randterme erreichen.

k_{ij} := k_{ij} + A_{Γ_{1}}(Φ_{i}, Φ_{j}) i, j = 1, …, M,
f_{i} := f_{i} + F_{Γ_{1}}(Φ_{i}) i = 1, …, M,

wobei K = (k_{ij})_{i,j=1,...,M}, F := (f_{i})_{i=1,...,M}. Wie schon zuvor kann hier elementweise vorgegangen werden. Es sollte beachtet werden, dass A_{Γ_{1}}(Φ_{i}, Φ_{j}) ≠ 0 nur gilt, wenn der i−te und der j−te Knoten zu einem gemeinsamen Elementrand auf Γ_{1} gehören. Ähnlich muss man bei F_{Γ_{1}}(Φ_{i}) nur Knoten auf Γ_{1} betrachten.

Zur Implementierung der Dirichlet-Randwerte gehen wir ebenfalls wie in der Variationsformulierung (1.4) vor. Wir konstruieren eine Funktion u_{g,h}, die u_{g} ∈ H_{Γ_{0},g_{0}}^{1}(Ω) möglichst gut approximiert. Dabei werden wir jedoch die Bedingung der Konformität verletzen, d.h. u_{g,h} ∉ H_{Γ_{0},g_{0}}^{1}(Ω). In der angelsächsischen Literatur wird ein solches Vorgehen nach Strang als ”varitional crime“ bezeichnet, vgl. Braess [2], Brenner/Scott [3]. Die Konstruktion wird dabei wie folgt vorgenommen

u_{g,h}(x, y) := ∑_{i∈J_{0}} g_{0}(x_{i}, y_{i})Φ_{i}(x, y) ∀(x, y) ∈ \overline{Ω}.

Offenbar gilt u_{g,h} ∈ H^{1}(Ω) und u_{g,h}(x_{i}, y_{i}) = g_{0}(x_{i}, y_{i}), i ∈ J_{0}. Allerdings interpoliert u_{g,h} nicht die Funktion g_{0} auf dem Randstück Γ_{0}. Somit gilt u_{g,h} ∉ H_{Γ_{0},g_{0}}^{1}(Ω). Es ist allerdings zu erwarten, dass u_{g,h} eine gute Approximation an u_{g} darstellt, sofern nur das Gitter fein genug gewählt wird.

In der Implementierung setzt man einfach die zu Randknoten gehörenden Einträge auf der Hauptdiagonalen der Steifigkeitsmatrix auf Eins und alle anderen Einträge, die zu dieser Spalte und Zeile gehören auf Null. Wie in (1.4) werden die inhomogenen Dirichlet–Randwerte dann in den Lastvektor integriert. Für alle j ∈ J_{0}, berechne

f_{i} := f_{i} − g_{0}(x_{i}, y_{i})k_{ij} ∀i ∈ J
k_{ij} := k_{ji} := 0 ∀i ∈ J
k_{ji} := 0 ∀i ∈ J_{0}, i ≠ j
f_{j} := g_{0}(x_{j}, y_{j})
k_{jj} := 1

Hiermit haben wir die komplette Struktur zur Programmierung der Methode der Finiten Elemente zur Verfügung. Was uns jetzt noch fehlt sind Routinen zur Berechnung der Elementmatrizen und –vektoren, sowie der Randintegrale zur Implementierung der inhomogenen Neumann–Randbedingungen.

Bei der Berechnung der Elementmatrizen/-vektoren werden wir ausnutzen, dass jedes Element affines Bild eines Referenzelements ist. In einem Element E gilt u_{h} = Φ^{t}u_{E}, wobei

Φ := [ Φ_{1} ; ⋮ ; Φ_{m} ] , u_{E} := [ u_{E,1} ; ⋮ ; u_{E,m} ] .

Daraus folgt ∇u_{h} = (∇Φ)^{t}u_{E} mit

∇Φ := [ ⋮  ⋮_{\frac{∂Φ_{m}}{∂x} \frac{∂Φ_{m}}{∂y}}^{\frac{∂Φ_{1}}{∂x} \frac{∂Φ_{1}}{∂y}} ] = [ (∇Φ_{1})^{t} ; ⋮ ; (∇Φ_{m})^{t} ] = [\frac{∂Φ}{∂x}, \frac{∂Φ}{∂y}] .

Man beachte, dass ∇Φ eine Matrix und kein Vektor ist. Mit diesen Bezeichungen ergibt sich

A_{E}(u_{h}, v_{h}) = ∫_{E} a∇u_{h}∇v_{h} + bu_{h}v_{h}dx
= u_{E}^{t} (∫_{E} a∇Φ(∇Φ)^{t}dx) v_{E} + u_{E}^{t} (∫_{E} bΦΦ^{t}dx) v_{E}

Offensichtlich gilt für die Elementmatrix

K_{E} = (∫_{E} a∇Φ(∇Φ)^{t} + bΦΦ^{t}dx) .

Für den Elementvektor ergibt sich analog

F_{E} = ∫_{E} fΦdx,

wobei f := [f_{1}, …, f_{m}] den Vektor darstellt, der die Werte der Funktion f_{h} in den Knoten darstellt.

Um die Berechnungen auf möglichst wenige Standardoperationen zurückzuführen, werden wir jetzt betrachten wie die Integration über ein beliebiges Element E durch Integration über ein Referenzelement Ê berechnet werden kann.

Sei nun E das Bild eines Referenzelements unter einer bijektiven, affinen linearen Abbildung T : Ê → E. Die Funktionen der nodalen Basis auf Ê und E haben folgende Beziehung zueinander:

Φ(x, y) = Φ̂(x̂, ŷ), mit (x, y) = T (x̂, ŷ) = (T_{1}(x̂, ŷ), T_{2}(x̂, ŷ)).

Für die Ableitungen der lokalen Basisfunktionen folgt mit der Kettenregel

\frac{∂Φ(x, y)}{∂x} = \frac{∂Φ(x, y)}{∂x̂} \frac{∂x̂}{∂x} + \frac{∂Φ(x, y)}{∂ŷ} \frac{∂ŷ}{∂x}
= \frac{∂Φ̂(x̂, ŷ)}{∂x̂} \frac{∂(T^{−1})_{1}(x, y)}{∂x} + \frac{∂Φ̂(x̂, ŷ)}{∂ŷ} \frac{∂(T^{−1})_{2}(x, y)}{∂x},

oder äquivalent

∇Φ(x, y) = ∇Φ̂(x̂, ŷ)(∇(T^{−1}))^{t}(x, y) = ∇Φ̂(x̂, ŷ)(∇T)^{−t}(x, y).

Für die Berechnung der Elementmatrizen und -vektoren ergibt sich somit

K_{E} = ∫_{E} ∇Φ̂K̂(∇Φ̂)^{t}dx̂dŷ + ∫_{E} b(T (x̂, ŷ))Φ̂Φ̂^{t}|det(∇T (x̂, ŷ))|dx̂dŷ,

wobei

K̂(x̂, ŷ) := a(T (x̂, ŷ))((∇T (x̂, ŷ))^{−t}(∇T (x̂, ŷ))^{−1})|det(∇T (x̂, ŷ))|.

Für den Elementvektor ergibt sich

F_{E} = ∫_{E} f(T (x̂, ŷ))Φ̂|det(∇T (x̂, ŷ)|dx̂dŷ.

Beispiel 1.3 Bei einer affin linearen Transformation T : Ê → E, T [x̂, ŷ]^{t} = B[x̂, ŷ]^{t}+ b gilt

∇T = B.

Um die Elementmatrizen und -vektoren auch im Fall variabler Koeffizienten a(·), b(·) berechnen zu können, benötigen wir noch passende Quadraturformeln. Hier bieten sich die Gauß’schen Integrationsformeln für Dreiecke und Rechtecke an, da sie für Polynome des jeweiligen Grades exakt sind. Im allgemeinen haben diese Formeln die Form

∫_{E} f(x̂, ŷ)dx̂dŷ ≈ ∑_{r=1}^{q} w_{r}f(b̂_{r}).

Dabei sind die Gewichte mit w_{1}, …, w_{r} und die Integrationspunkte mit b̂_{1}, …, b̂_{r} bezeichnet.

Beispiel 1.4 Für linear Dreiecke (Lagrange) mit dem Einheitsdreieck als Referenzelement Ê ist der Integrationspunkt b̂_{1} = b̂ := (1/3, 1/3) und das Gewicht ist w_{1} = w = |Ê| = 1/2. Die so erhaltene Quadraturformel ist exakt für lineare Funktionen.

Wir geben nun eine algorithmische Beschreibung zur Aufstellung der Steifigkeitsmatrix K aus den Elementmatrizen an. Dabei bezeichnen wir mit k_{rs}^{(l)} die Koeffizienten der Elementsteifigkeitsmatrix K_{E_{l}}, d.h. K_{E_{l}} =: (k_{rs}^{(l)})_{r,s=1,...,m}.

Algorithmus 1.1 Berechnung der Steifigkeitsmatrix K ohne Berücksichtigung der Randbedingungen Initialisierung: K = O Berechnung K: for l := 1 : L

for r := 1 : m for s := 1 : m Berechne k_{rs}^{(l)} i := i_{r}^{(l)} := ELNODE[l,r] j := j_{s}^{(l)} := ELNODE[l,s] k[i, j] := k[i, j] + k_{rs}^{(l)}

endfor endfor endfor

Zur Einarbeitung der Dirichlet–Randbedingungen muss man nur die Dirichlet– Randknoten auf Γ_{0} durchlaufen, d.h. die Indexmenge J_{0}, und wie zuvor beschrieben vorgehen. Die dazu benötigte Liste BONODE wird während der Gittererzeugung angelegt, s. Abschnitt 1.8.1.

Schließlich müssen wir noch die inhomogenen Neumann– bzw. Robin– Randbedingungen programmieren. Der Einfachheit halber beschränken wir uns dazu auf lineare Dreieckselemente. Bei anderen Elementen kann jedoch prinzipiell genauso verfahren werden.

Wir betrachten zunächst

A_{Γ_{1}}(Φ_{i}, Φ_{j}) = ∫_{Γ_{1}} a c Φ_{i}Φ_{j}dσ.

Zur Berechnung dieses Integrals müssen wieder alle Knoten auf Γ_{1} besucht werden. Dieses geschieht wie zuvor elementweise. In der Liste EDGE (vgl. Abschnitt 1.8.1), findet man alle Elemente, die auf dem Rand Γ_{1} liegen, sowie deren Anfangsknoten. Aus dieser Information lassen sich mit Hilfe der lokalen Knotennummerierung die weiteren Knoten auf dem Rand des Elements bestimmen. In unserem Fall, dem linearer Dreieckselemente, ist es nur ein weiterer Knoten.

Wir betrachten nun ein beliebiges Randelement E_{l} ⊂ Γ_{1} und der Anfangsknoten sei mit r_{1}, der weitere Knoten mit r_{2} bezeichnet. Benutzen wir wieder die Relation zwischen lokaler und globaler Knotennummerierung, die durch e_{i}^{l} gegeben ist, so können wir die Randknoten wie folgt darstellen

(x_{i_{1}}, y_{i_{1}}), (x_{i_{2}}, y_{i_{2}}) mit i_{1} := e_{r_{1}}^{l}, i_{2} := e_{r_{2}}^{l}.

Wir benötigen außerdem die Indexmengen

Ŝ_{t} := (r_{1}, r_{2}), S_{t} := {(r, s) : r, s ∈ Ŝ_{t}, r ≥ s}.

Bei linearen Dreiecken kann A_{Γ_{1}}(Φ_{i}, Φ_{j}) durch Zusammensetzen zweier Integrale der Form

I_{r,s}^{(t)} := ∫_{(x_{i_{1}},y_{i_{1}})}^{(x_{i_{2}},y_{i_{2}})} a c Φ_{i}Φ_{j}dσ, (r, s) ∈ S_{t},

erfolgen.

Wie schon bei der Berechnung der Steifigkeitsmatrix ohne (inhomogene) Randbedingungen, kann man nun elementweise vorgehen, etwa wie folgt

∀t ∈ {1, …, L_{1}}
∀(r, s) ∈ S_{t}
Berechne I_{rs}^{(t)}
i := e_{r}^{t}, j := e_{s}^{t}
k_{ij} := k_{ij} + I_{rs}^{(t)}, k_{ji} := k_{ij}, ∀i ≠ j
k_{ii} := k_{ii} + I_{rs}^{(t)}, ∀i = j

Um I_{rs}^{(t)} zu berechnen, wenden wir wie zuvor die affin-lineare Transformation T vom Einheitsdreieck auf das jeweilige Element an. Die Randintegrale müssen dann nur noch auf dem Rand des Einheitsdreiecks durchgeführt werden. Dabei treten drei verschiedene Fälle auf. Bezeichnen wir mit

d := ((x_{i_{2}} − x_{i_{1}})^{2} + (y_{i_{2}} − y_{i_{1}})^{2})^{1/2}

die Länge des Randes im Originalelement und führen zur Abkürzung die Funktion

Ĝ(x̂, ŷ) := a(T (x̂, ŷ))c(T (x̂, ŷ))

ein, so gilt es die drei folgenden verschiedenen Grundintegrale zu berechnen

(i) d ∫_{0}^{1} Ĝ(x̂, 0)Φ̂_{r}(x̂, 0)Φ̂_{s}(x̂, 0)dx̂,

(ii) \frac{d}{\sqrt{2}} ∫_{0}^{√\overline{2}} Ĝ(x̂, ŷ)Φ̂_{r}(x̂, ŷ)Φ̂_{s}(x̂, ŷ)dσ̂, x̂ = −\frac{ŝ}{\sqrt{2}} + 1, ŷ = \frac{ŝ}{\sqrt{2}},

(iii) d ∫_{0}^{1} Ĝ(0, ŷ)Φ̂_{r}(0, ŷ)Φ̂_{s}(0, ŷ)dŷ.

Zur numerischen Integration benutzen wir eindimensionale Gauß’sche Quadraturformeln mit den Integrationspunkten b̂_{1}, b̂_{2} und dem Gewicht w_{1}:

∫_{a}^{b} f(η)dη ≈ w_{1} [f(b̂_{1}) + f(b̂_{2})] ,

wobei

b̂_{1} := −\frac{1}{2\sqrt{3}}(b − a) + \frac{b + a}{2}, b̂_{2} := \frac{1}{2\sqrt{3}}(b − a) + \frac{b + a}{2}, w_{1} := \frac{b − a}{2}.

Es bleibt jetzt nur noch

F_{Γ_{1}}(Φ_{i}) = ∫_{Γ_{1}} a g_{1} Φ_{i}dσ, ∀i = 1, …, M,

zu berechnen. Dieses Integrale kann wieder durch Zusammensetzen von Integralen der Form

I_{r}^{t} := ∫_{(x_{i_{1}},y_{i_{1}})}^{(x_{i_{2}},y_{i_{2}})} a g_{1} Φ_{t}dσ, r ∈ Ŝ_{t}.

Dabei können wir wie oben beschrieben vorgehen.

∀t ∈ {1, …, L_{1}}
∀r ∈ Ŝ_{t}
Berechne I_{r}^{(t)}
i := e_{r}^{t}
f_{i} := f_{i} + I_{r}^{(t)}

Die Berechnung des Integrals I_{r}^{t} mit Hilfe der Transformation auf das Referenzdreieck und unter Benutzung numerischer Ingration kann wie oben angegeben durchgeführt werden.

18

Literaturverzeichnis

[1] O. Axelsson and V.A. Barker. Finite Element Solution of Boundary Value Problems: Theory and Computation. Academic Press, 1984.

[2] Dietrich Braess. Finite Elemente. Springer-Verlag, 2nd edition, 1997.

[3] Susanne C. Brenner and L. Ridgway Scott. The Mathematical Theory of Finite Element Methods. Springer-Verlag, 1994.

[4] Philippe G. Ciarlet. The Finite Element Method for Elliptic Problems. North-Holland, 1978.

[5] C.W. Cryer. Numerik Partieller Differentialgleichungen II. Vorlesungsskript Sommersemester 1998.

[6] P.L. George. Automatic Mesh Generation: Application to Finite Element Methods. John Wiley & Sons, 1991.

[7] W. Hackbusch. Theorie und Numerik elliptischer Differentialgleichungen. Teubner-Verlag, Stuttgart, 1986. Engl. Transl.: Elliptic Differential Equations. Theory and Numerical Treatment. Springer, New York, 1992.

[8] J.L. Lions and E. Magenes. Non–Homogeneous Boundary Value Problems and Applications, Part I. Springer–Verlag Berlin Heidelberg New York, 1972.

[9] J. Nečas. Les Méthodes Directes en Théorie des Equations Elliptiques. Masson-Paris/Academia-Prague, 1967.

[10] J. Wloka. Partielle Differentialgleichungen. B.G. Teubner Stuttgart, 1982.