4.6.1 Allgemeine Vorgehensweise
Zunächst muss in geeigneter Form das Problem beschrieben werden, nämlich die Differentialgleichung, die Randbedingungen und das Gebiet. Ebenfalls zum Preprocessing gehört die Vernetzung des Gebiets. Je allgemeiner die zugelassenen Gebiete sind, desto höher sind die Anforderungen an den Vernetzungsalgorithmus.
Auch wenn man bei der Aufstellung des Gleichungssystems an knotenweise definierte Funktionen denkt, wird die Steifigkeitsmatrix und die rechte Seite elementweise erzeugt und assembliert. Die Integration wird dabei näherungsweise mit Quadraturformeln ausgeführt. Der dabei entstehende Fehler kann mit Hilfe des 2. Lemmas von Strang abgeschätzt werden, wird jedoch in den theoretischen Betrachtungen meist vernachlässigt.
Ein Vorteil der Finite-Elemente-Methode besteht darin, dass die Steifigkeitsmatrix nur dünn besetzt ist (sparse). Für deren Speicherung stehen verschiedene Techniken zu Verfügung. Verbreitet ist die Kompaktliste zeilenweise (compressed row storage), wobei nur die Nichtnullelemente der Matrix gespeichert werden.
Das entstehende Gleichungssystem kann mit einem direkten Verfahren gelöst werden, das die dünne Besetztheit geschickt ausnutzt. Bei großer Dimension des Gleichungssystems sind jedoch iterative Verfahren wie das Verfahren der konjugierten Gradienten, effektiver, vor allem wenn sie mit einer guten Vorkonditionierung kombiniert werden. Die besten Vorkonditionierungen basieren auf einer Netzhierarchie, die die zu implementierende Software entsprechend komplex werden lässt.
Schließlich gehört zu jeder Implementierung ein Postprocessing, das die graphische Darstellung der Ergebnisse, die Berechnung mittelbarer Größen (Ableitungen oder Funktionalwerte der Lösung) und eine a posteriori Fehlerbewertung des Ergebnisses beinhaltet. An dieser Stelle entsteht im Allgemeinen eine Schleife, wobei man aus den a posteriori Fehlerinformationen einen neuen, besser geeigneten Finite-Elemente-Raum konstruiert und wieder in die Phase der Assemblierung und Lösung des Gleichungssystems eintritt.
Auf Grund der Spezifik unterschiedlicher Typen partieller Differentialgleichungen und der Vielfalt möglicher Diskretisierungen und anderer Einzelbausteine (Speichertechniken, Löser, etc.) gibt es keine universellen Softwarelösungen.
4.6.2 Datenstrukturen
Ein Finite-Elemente-Netz besteht im 2D-Fall aus Knoten (vertices), Kanten
(edges) und Elementen
. Im 3D-Fall kommen noch Seiten
(faces) hinzu.
In jedem Fall gibt es ein Feld , in dem zu jedem Knoten die Koordinaten gespeichert sind.
Im einfachen Fall (reduzierte Datenstruktur) gibt es nur noch ein Feld , in dem zu jedem Element die Knotennummern gespeichert sind.
Im allgemeinen Fall (vollständige Datenstruktur) gibt es folgende Felder
2D:
3D:
Weitere Felder können nützlich sein:
- Liste der Randseiten
- Liste der Nachbarelemente
- Vater-Sohn-Beziehungen in hierarchischen Netzen
- Materialbeschreibung
Beispiel: Dreiecksnetz für
Listen:
entspricht einer Zuordnung von lokalen Knotennummern zu globalen Knotennummern. Diese ist nicht eindeutig.
4.6.3 Aufbau des Finite-Elemente-Gleichungssystems
Zu lösen ist das lineare Gleichungssystem mit
Die Berechnung der Bilinearform wurde bereits am Beispiel erläutert. Obwohl hier knotenbezogene Größen eingehen (), läuft die Schleife über die Elemente:
Für ein Element ist das Integral nur ungleich 0, wenn
und
Knoten von
sind.
Wir betrachten das lineare Dreieckselement mit den Ecken , und . Der Beitrag von zur Steifigkeitsmatrix hat die Form
(überall anders stehen Nullen!) Die auf Nicht-Nullelemente reduzierte Matrix
heißt Elementsteifigkeitsmatrix. Die Beziehung zwischen und
wird durch
definiert. Die Addition der Elemente von
an die entsprechenden Stellen von
heißt Assemblierung.
Auch für den Lastvektor (rechte Seite) führen wir die Berechnung am Referenzelement durch und transformieren dann auf das echte Element:
Die Daten werden also nur an Punkten ausgewertet. Dies wird bei stark oszillierenden Funktionen kritisch (Shannon-Theorem!).
4.6.4 Berechnung der Element-Steifigkeitsmatrix
Die Integrale berechnet man mit Hilfe von Quadraturformeln, die Quadraturpunkte über das Referenzsystem:
Einbau von homogenen Dirichlet-Randbedingungen: Man assembliert die Matrix und die rechte Seite und streicht dann Zeilen und Spalten, die zu Randknoten gehören.
Eigenschaften der Systemmatrix:
Die Matrix
- hat die Dimension
- ist dünn besetzt, die Anzahl der Nicht-Nullelemente pro Zeile ist
, die der gesamten Matrix
- ist schlecht konditioniert:
unabhängig von der Raumdimension
4.6.5 Methoden zur Netzgenerierung
4.6.5.1 Manuelle Methode
Eingabe aller Daten per Hand. Das ist praktikabel für Gebiete mit einfacher Geometrie und sehr große Netze.
4.6.5.2 Produktmethode
Geeignet für prismatische oder zylinderförmige Gebiete:
4.6.5.3 Transformation von Referenznetzen
Nützlich bei topologisch einfachen Gebieten mit gekrümmten Rändern.
4.6.5.4 Halbautomatische Methode
Vorgehen:
- Zerlegung des Gebiets in Teilgebiete
- Teilgebiete werden mit einer der vorherigen Methoden vernetzt
Achtung: Sichern einer zulässigen Vernetzung!
4.6.5.5 Quadtree- und Octree-Methode
Vorgehen:
- Start mit 1 Rechteck, das
enthält
- rekursive Teilung der Elemente, bis der Rand von
hinreichend gut approximiert wird
- Streichen aller Elemente, die außerhalb von
liegen
- Verschieben von Knoten und Kanten zur besseren Randapproximation
Man kann hier sehr schnell bestimmen, in welchem Element ein gegebener Punkt liegt.
4.6.5.6 Delauney-Vernetzung
Gegeben: Punkte
Vorgehen:
- Konstruiere Voronoi-Polyeder
durch Mittelsenkrechten
- Verbinde diejenigen Punkte, deren Voronoi-Polyeder eine gemeinsame Kante oder Fläche besitzen
Die Delauney-Triangulierung trianguliert die konvexe Hülle der Punkte .
Im Inneren jeden Umkreises zu einem Dreieck liegt kein weiterer Punkt. D.h. die Summe der einer Seite gegenüberliegenden Innenwinkel ist kleiner als oder gleich 180°.
4.6.5.7 Advancing-Front-Methoden
Vorgehen:
- Vernetze den Gebietsrand
- Generiere ein Element, das mindestens eine Seite auf der Front besitzt, nach verschiedenen Kriterien, z.B.
- Lücken schließen
- spitze Winkel in Front eliminieren
- Definiere neue Front (Rand des noch nicht vernetzten Gebietes). ist diese leer, gehe zum zweiten Schritt.
4.6.6 Netzverfeinerung
Man kann geeignete Finite-Elemente-Netze auch durch sukzessive Verfeinerung von groben Ausgangsnetzen gewinnen. Die lokale Netzdichte kann gesteuert werden durch
- a-priori Informationen
- a-posteriori Informationen (lokale Fehlerindikatoren, die basierend auf der Lösung auf den bisherigen Netz berechnet werden)
Teilung von Elementen:
2D: Viertelung von Dreiecken (rote Verfeinerung) erzeugt kongruente Dreiecke
3D: Achtelung von Tetraedern liefert im Allgemeinen keine kongruenten Tetraeder. Es existieren verschiedene Strategien.
Falls nicht alle Dreiecke geviertelt werden, entsteht im Allgemeinen keine zulässige Zerlegung (hängende Knoten entstehen). Abhilfe: angrenzende Dreiecke werden halbiert (grüne Verfeinerung)
Um nicht immer kleinere Winkel zu erhalten, werden grüne Dreiecke nicht nochmals halbiert.