\noindent
Sei nun $\T_{\ell} = \{T_1,T_2,\ldots,T_N\}$ eine Partition des Randes $\Gamma$ in $N$ Randstücke mit den charakteristischen Funktionen
-\begin{align}
+\begin{align*}
\chi_i(\bs x) &=
\begin{cases}
1&\text{für alle }\bs x \in T_i\\
0&\text{sonst}
\end{cases}
\quad \text{für alle } i\in\{1,2,\ldots,N\},
-\end{align}
-und sei $\P^0(\T_{\ell}) := span\{\chi_i : i=1,\ldots,N\}$ der aufgespannte Teilraum von $\tilde H^{-1/2}(\Gamma)$. Wir suchen also die Lösung $\phi_{\ell} \in \P^0(\T_{\ell})$ für
-\begin{align}
- \llangle \phi_{\ell},\psi_{\ell} \rrangle & = \langle f,\psi_{\ell} \rangle\quad\text{für alle } \psi_{\ell} \in \P^0(\T_{\ell}),
-\end{align}
+\end{align*}
+und sei $P^0(\T_{\ell}) := span\{\chi_i : i=1,\ldots,N\}$ der aufgespannte Teilraum von $\tilde H^{-1/2}(\Gamma)$. Wir suchen also die Lösung $\phi_{\ell} \in P^0(\T_{\ell})$ für
+\begin{align*}
+ \llangle \phi_{\ell},\psi_{\ell} \rrangle & = \langle f,\psi_{\ell} \rangle\quad\text{für alle } \psi_{\ell} \in P^0(\T_{\ell}),
+\end{align*}
was aufgrund der Basiseigenschaft von $\chi_i$ eindeutig und äquivalent ist zu
-\begin{align}
+\begin{align}\label{math:slp:gls:galerkin}
\llangle \phi_{\ell},\chi_j \rrangle & = \langle f,\chi_j \rangle\quad\text{für alle } j\in\{1,2,\ldots,N\}.
\end{align}
Dadurch können wir das Problem als Summe anschreiben
-\begin{align}
+\begin{align*}
\sum_{j=1}^N x_{\ell,j} \llangle \chi_j,\chi_i\rrangle = \langle g, \chi_i\rangle,
-\end{align}
+\end{align*}
wobei $\bs x_{\ell} = \{x_{\ell,1},x_{\ell,2},\ldots,x_{\ell,N}\}$ der Koordinatenvektor von $\phi_{\ell}$ zur Basis $\chi_1,\ldots,\chi_N$ ist.
So schreiben wir
\begin{align}
\begin{align*}
L_j(x) &= \prod_{i=0 \atop i\neq j}^p \frac{x-x_i}{x_j-x_i}.
\end{align*}
-Mit diesem Wissen definieren wir uns nun den Interpolationsoperator $\I_p : \C[0,1]\to\P^{p}$
+Mit diesem Wissen definieren wir uns nun den Interpolationsoperator $\I_p : \C[0,1]\to P^{p}$
\begin{align*}
\I_pu &:= \sum_{j=0}^p u(x_j)L_j.
\end{align*}
A_{jk} = \int_{T_j} \int_{T_k} \frac{1}{|\bs x- \bs y|} ds_{\bs y} ds_{\bs x} \in \R^3.
\end{align}
auf zwei beschränkten, achsenorientierten Rechtecken $T_j,T_k \subseteq\R^3$ beschäftigen. Die im Folgenden auftretenden Stammfunktionen $\int f(x) dx$ werden wir der Einfachheit halber jeweils mit additiver Verschiebung $0$ schreiben.
-Dazu wollen wir \cite{mai:3dbem} folgend, zwei Stammfunktionen zitieren, welche durch das Aufspalten des Integrals \eqref{math:analy:int} auftreten werden.
+Dazu wollen wir, \cite{mai:3dbem} folgend, zwei Stammfunktionen zitieren, welche durch das Aufspalten des Integrals \eqref{math:analy:int} auftreten werden.
\\\noindent
\begin{lem}
Für Die Stammfunktion
\end{align*}
benötigen, welche alle in \cite[Seite 2-3]{mai:3dbem} bewiesen wurden.
\end{lem}
-Hierbei wollen wir darauf hinweisen, dass wir für $p=-1$ nicht die in \cite[Seite 3]{mai:3dbem} vorgeschlagene Formel verwenden, sondern die selbst hergeleitete. Mithilfe von Substitution durch $z = \frac{y-x}{\abs{\lambda}}$ und $dz = \frac 1 {\abs{\lambda}} dy$ gilt
+Hierbei wollen wir darauf hinweisen, dass wir für $p=-1$ nicht die in \cite[Seite 3]{mai:3dbem} vorgeschlagene Formel verwenden, sondern eine selbst hergeleitete. Mithilfe von Substitution durch $z = \frac{y-x}{\abs{\lambda}}$ und $dz = \frac 1 {\abs{\lambda}} dy$ gilt
\begin{align*}
\int \frac 1 {(x-y)^2 +\lambda^2} dy &= \int \frac 1 {\lambda^2((\frac{y-x}{\abs{\lambda}})^2 +1)} dy
= \frac 1 {\abs{\lambda}} \int \frac 1 {z^2+1} dz\\
\end{align*} welche ebenfalls in \cite[Seite 5-7]{mai:3dbem} bewiesen wurden.
\end{lem}
-\subsection{Integral über zwei Elemente}
-Bei der Berechnung von \eqref{math:analy:int} werden wir Aufgrund der speziellen Form des Kerns geometrisch zwischen zwei Fällen zu unterscheiden. Entweder die beiden Elemente liegen geometrisch in parallelen Ebenen oder in orthogonalen Ebenen.
+\subsection{Integral über zwei Elemente}\label{sec:analyt:int}
+Bei der Berechnung von \eqref{math:analy:int} werden wir Aufgrund der speziellen Form des Kerns geometrisch zwischen zwei Fällen unterscheiden. Entweder die beiden Elemente liegen geometrisch in parallelen Ebenen oder in orthogonalen Ebenen.
% \subsubsection{Parallele Elemente}
\noindent
\begin{align*}
T_j &= \{\bs v + \lambda_1 a {\bs a} + \lambda_2 b {\bs b} ~|~ \lambda_1,\lambda_2 \in[0,1]\}\\
T_k &= \{\tilde {\bs v} + \tilde \lambda_2 \tilde b \tilde {\bs b} + \tilde \lambda_1 \tilde a \tilde {\bs a} ~|~ \tilde \lambda_1,\tilde \lambda_2 \in[0,1]\}
-\end{align*} wobei wir hier ${\bs b} = \tilde {\bs a}$ und ${\bs b} \neq \tilde{\bs a}$ annehmen. Wir setzen $\bs{\delta} = \tilde \v - \v \in \R^3$.
+\end{align*} wobei wir hier ${\bs b} = \tilde {\bs a}$ und ${\bs a} \perp \tilde{\bs b}$ annehmen. Wir setzen $\bs{\delta} = \tilde \v - \v \in \R^3$.
Dann gilt nach \cite{mai:3dbem}
\begin{align*}
&\int_{T_j} \int_{T_k} \frac{1}{|\bs x- \bs y|} ds_{\bs y} ds_{\bs x}\\
\subsection{Bestimmtes Integral}
-Mithilfe der Stammfunktionen $F_p$ und $F_o$ können wir nun das Integral aus \eqref{math:analy:int} für orthogonal und parrallel liegenden achsenorientierte Rechtecke $T$ und $\tilde T$ berechnen.
+Mithilfe der Stammfunktionen $F_p$ und $F_o$ können wir nun das Integral aus \eqref{math:analy:int} für orthogonal und parrallel liegenden achsenorientierte Rechtecke $T$ und $\tilde T$ berechnen. Im Folgenden sei $F_{p/o} \in \{F_p,F_o\}$.
\begin{align*}
\int_T& \int_{\tilde T} \frac{1}{|x-y|} ds_y ds_x\\
&= \int_0^{k_1}\int_0^{k_2}\int_0^{\tilde k_1}\int_0^{\tilde k_2}
Weiterhin gehen wir auch auf die Implementierung des Verfeinern Algorithmus aus Kapitel \ref{sec:bem:ref} ein.
Dabei werden wir wichtige Schritte anhand des Codes und kleinen Beispielen hervorheben.
-\subsection{Datenstruktur}
+\subsection{Datenstruktur}\label{sec:implement:daten}
Für die Implementierung in \Matlab~und \Cpp~werden wir eine einheitliche Datenstruktur einführen.
Die für die Partition $\T_{\ell} = \{T_1\ldots T_N\}$ benötigten Knoten $\K_{\ell} = \{k_1\ldots k_M\}$ stellen wir in einer $M \times 3$ Matrix $COO$ dar. Dabei enthält die $j$-te Zeile die Koordinaten des Knoten $C_j$, d.h. :
\begin{align}
\begin{align}
ELE[i,1:4] = T_i := (j,k,l,m).
\end{align}
-Die Knoten ordnen wir gegen den Uhrzeigersinn an und der Knoten $C_j$ sei der Punkt $\bs v$ aus Definition \ref{thm:def:T}, also der kleinste Bezüglich des Koordinatensystems.
+Die Knoten ordnen wir gegen den Uhrzeigersinn an und der Knoten $C_j$ sei der Punkt $\bs v$ aus Definition \ref{thm:def:T}.
\noindent
-Für die bessere Handhabung der Elemente beim Verfeinern der Partition, müssen wir auch die Nachbarschaftsrelationen geeignet abspeichern. Aufgrund der Netzstabilität werden wir maximal zwei Nachbarn pro Kante eines Elements zulassen.
+Für die bessere Handhabung der Elemente beim Verfeinern der Partition, müssen wir auch die Nachbarschaftsrelationen geeignet abspeichern.
+Im Folgenden werden wir maximal zwei Nachbarn pro Kante eines Elements zulassen, damit das Verhältnis der Elementgrößen benachbarter Elemente beschränkt bleibt.
Wir legen also eine $M \times 8$ Matrix für die Indizes der Nachbarelemente an, wobei die $i$-te Zeile die Nachbarelemente $\{T_{n_1},\ldots,T_{n_8}\}$ zum Element $T_i$ enthält.
\begin{align}
NEI[i,1:8] = N_i := (n_1,\ldots,n_8)
\end{align}
-Offensichtlich ist $i \notin N_i$. Wir wollen uns aber noch genauer eine geeignete Anordnung für die Nachbarelemente überlegen. Hierbei bezeichnen wir die Seite $[j,k]$ eines Elements als Seite 1 und gegen den Uhrzeigersinn alle weiteren mit $\{2,3,4\}$. Für den einfacheren Zugriff auf die Elemente einer Seite $s$, seien die Nachbarelemente zur Seite $s$ unter den Indizes $n_s$ und $n_{s+4}$ abgespeichert. Die Nachbarelemente $\{T_{n_s}, T_{n_{s+4}}\}$ liegen also an der Seite $s$ des Elements. Für Seiten die nur einen Nachbarn $T_{n_s}$ besitzen setzen wir $n_{s+4}=0$ und für Seiten mit keinem Nachbarn setzen wir $n_s = n_{s+4} = 0$. Daraus folgt unmittelbar, dass für $n_s =0$ auch $n_{s+4} = 0$ gilt, die Seite $s$ also keine Nachbarelemente besitzt und umgekehrt folgt aus $n_{s+4} \neq 0$ $n_s \neq 0$, womit die Seite $s$ genau zwei Nachbarelemente hat.
-(Siehe Abb.:\ref{exmpl13:nei:part})
+Offensichtlich ist $i \notin N_i$. Wir wollen uns aber noch genauer eine geeignete Anordnung für die Nachbarelemente überlegen. Hierbei bezeichnen wir die Seite $[j,k]$ eines Elements als Seite 1 und gegen den Uhrzeigersinn alle weiteren $[k,l], [l,m]$ und $[m,j]$ mit 2,3,4. Für den einfacheren Zugriff auf die Elemente einer Seite $s \in \{1,2,3,4\}$, seien die Nachbarelemente zur Seite $s$ unter den Indizes $n_s$ und $n_{s+4}$ abgespeichert. Die Nachbarelemente $\{T_{n_s}, T_{n_{s+4}}\}$ liegen also an der Seite $s$ des Elements. Für Seiten die nur einen Nachbarn $T_{n_s}$ besitzen setzen wir $n_{s+4}=0$ und für Seiten mit keinem Nachbarn setzen wir $n_s = n_{s+4} = 0$. Daraus folgt unmittelbar, dass für $n_s =0$ auch $n_{s+4} = 0$ gilt, die Seite $s$ also keine Nachbarelemente besitzt und umgekehrt folgt aus $n_{s+4} \neq 0$ $n_s \neq 0$, womit die Seite $s$ genau zwei Nachbarelemente hat.
+(Siehe Abbildung \ref{exmpl13:nei:part})
\begin{figure}[ht]
\centering
\subfloat[Lage]{\includegraphics[width=0.4\textwidth]{fig/exmpl13_nei_part}}
-\subfloat[Nachbarn]{\input{fig/exmpl13_nei_part}}
-\caption{\scriptsize An dieser Stellen wollen wir die Nachbarschaftsrelationen des Elements 4 aus Abb.\ref{exmpl13} hervorheben, denn das Element hat an den ersten beiden Seiten jeweils einen Nachbarn, an der dritten Seite zwei Nachbarn und an der vierten Seite keinen. Der Index der Nachbarelemente ist jeweils in den Nachbarschaftsrelationen gespeichert. Wobei bei Seiten mit nur einem Nachbar der Nachbar immer im kleineren Index gespeichert wird. Dass das Element 4 selbst Nachbar einer Seite mit zwei Elementen ist geht aus den Informationen für Element 4 nicht hervor, ist aber dafür in Element 9 vermerkt. Sollte beispielsweise das Element 7 oder 12 horizontal geteilt werden, so muss klarer Weise auch Element 4 mindestens horizontal geteilt werden, welches auch eine horizontale Teilung von Element 9 erzwingen würde.}
+\subfloat[Nachbarn : vierte Zeile der NEI Matrix]{\input{fig/exmpl13_nei_part}}
+\caption{\scriptsize An dieser Stellen wollen wir die Nachbarschaftsrelationen des Elements 4
+% aus Abb.\ref{exmpl13}
+hervorheben, denn das Element hat an den ersten beiden Seiten jeweils einen Nachbarn, an der dritten Seite zwei Nachbarn und an der vierten Seite keinen. Der Index der Nachbarelemente ist jeweils in den Nachbarschaftsrelationen gespeichert. Wobei bei Seiten mit nur einem Nachbar der Nachbar immer im kleineren Index gespeichert wird. Dass das Element 4 selbst Nachbar einer Seite mit zwei Elementen ist geht aus den Informationen für Element 4 nicht hervor, ist aber dafür in Element 9 vermerkt. Sollte beispielsweise das Element 7 oder 12 horizontal geteilt werden, so muss klarer Weise auch Element 4 mindestens horizontal geteilt werden, welches auch eine horizontale Teilung von Element 9 erzwingen würde.}
\label{exmpl13:nei:part}
\end{figure}
\begin{enumerate}
\item keine Teilung
\item volle Teilung in vier gleich große Elemente
-\item halbe Teilung in zwei gleichgroße vertikal getrennte Elemente
-\item halbe Teilung in zwei gleichgroße horizontal getrennte Elemente
+\item halbe Teilung in zwei gleich große vertikal getrennte Elemente
+\item halbe Teilung in zwei gleich große horizontal getrennte Elemente
\end{enumerate}
-Zusätzlich wurde auch Typ 5. belegt, welcher als Ergebnis eine volle Teilung 2. ausführt, diese aber schrittweise durch eine 3. Teilung und anschließend durch jeweils eine 4. Teilung. Aus Sicherheitsgründen wird auch jede vierte volle 2. Teilung durch eine 5. Teilung ausgeführt, da sonst kurzzeitig Seiten mit mehr als zwei Elementen auftreten könnten.
+Zusätzlich wurde auch Typ 5 belegt, welcher als Ergebnis eine volle Teilung vom Typ 2 ausführt, diese aber schrittweise durch eine Teilung vom Typ 3 und anschließend durch jeweils eine Typ 4 Teilung. Aus Sicherheitsgründen wird auch jede vierte Teilung vom Typ 2 durch eine Typ 5 Teilung ausgeführt, da sonst kurzzeitig Seiten mit mehr als zwei Elementen auftreten könnten.
\noindent
Damit jedem Element $T_i$ eine Teilungsart zugeordnet werden kann, haben wir einen Mar\-kierungs\-vektor $marked \in \{1,2,3,4,5\}^M$ eingeführt. Dabei entspricht $marked_i$ der Art der Teilung für das Element $T_i$. Um isotrope und auch uniforme Teilungen zu erleichtern kann statt dem Vektor $marked$ auch nur ein Skalar übergeben werden $marked \in {1,2\dots5}$, wodurch jedes Element mit der gewählten Art verfeinert wird.
Relevant zum Verfeinern eines Netzes sind also die Koordinaten $COO$, Elemente $ELE$ sowie die Nachbarschaftsrelationen $NEI$ und der Markierungs\-vektor $marked$.
\noindent
-Da wir später den Fehlerschätzer berechnen wollen, ist es wichtig sich zu jedem Element seine Teilelemente zu merken. Dazu legen wir während der Teilung eine $M \times 4$ Matrix an, in der die maximal vier Elementindizes gespeichert werden. Wenn wir also ein Element in vier gleich große Teile verfeinern, so wird das neue Element links unten das erste sein und alle weiteren folgen gegen den Uhrzeigersinn. Teilen wir ein Element in zwei gleich große Elemente, so werden die doppelt belegten Quadranten auch doppelt eingetragen. Ein gar nicht geteiltes Element wird also vier mal den alten Indizes speichern. Dadurch wird sichergestellt, dass das arithmetische Mittel über die Elemente immer gültig auszuführen ist.
+Da wir später den Fehlerschätzer berechnen wollen, ist es wichtig sich zu jedem Element seine Teilelemente zu merken. Dazu legen wir während der Teilung eine $M \times 4$ Matrix $F2S$ an, in der die maximal vier Elementindizes gespeichert werden. Wenn wir also ein Element in vier gleich große Teile verfeinern, so wird das neue Element links unten das erste sein und alle weiteren folgen gegen den Uhrzeigersinn. Teilen wir ein Element in zwei gleich große Elemente, so werden die doppelt belegten Quadranten auch doppelt eingetragen. Ein gar nicht geteiltes Element wird also vier mal den alten Indizes speichern. Dadurch wird sichergestellt, dass das arithmetische Mittel über die Elemente immer gültig auszuführen ist.
(Siehe Abbildung \ref{exmpl13})
\begin{lstlisting}[language=M, numbers=none]
[COO_fine, ELE_fine, NEI_fine, F2S] = refineQuad(COO, ELE, NEI, marked);
\label{exmpl1:f2s:part}
\end{figure}
-\subsection[Berechnung der A Matrix]{Berechnung der $A$ Matrix}
+\subsection[Berechnung der A Matrix]{Berechnung der $A$ Matrix}\label{sec:implement:A}
An dieser Stelle wenden wir uns noch ein mal der Berechnung der Matrix $A \in \R^{N\times N}$ zu, deren Einträge bestimmt sind durch das Integral
\begin{align} \label{math:imple:Ajk}
A_{jk} = \int_{T_j} \int_{T_k} \frac 1 {4\pi} \frac 1 {\abs{\bs x - \bs y}} ds_{\bs y} ds_{\bs x}.
\noindent
Innerhalb von \Matlab~kann der Code dann kompiliert werden
\begin{lstlisting}[language=M, numbers=none]
-mex -o 'mex_build_V.cpp' 'slpRectangle.cpp' CFLAGS="\$CFLAGS -fopenmp"...
- CXXFLAGS="\$CXXFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp"
+mex 'mex_build_V.cpp' 'slpRectangle.cpp' CFLAGS="\$CFLAGS -O2 -fopenmp"...
+ CXXFLAGS="\$CXXFLAGS -O2 -fopenmp" LDFLAGS="\$LDFLAGS -O2 -fopenmp"
\end{lstlisting}
Nach erfolgreichem Kompilieren wurde eine \lstinline!mex_build_V.mexa64! Datei für ein 64bit System erstellt, welche dann in \Matlab~über
\begin{lstlisting}[language=M, numbers=none]
A = mex_build_V(coo,ele,zeta,typ);
\end{lstlisting}
-aufgerufen werden kann. Hierbei steht wie in der Datenstruktur beschrieben \lstinline!coo! für die Koordinatenmatrix und \lstinline!ele! für die Elementmatrix. Weiterhin sei \lstinline!zeta! ein Array mit den Einträgen \lstinline!zeta! $=(q,\zeta_0,\zeta_1)$. Mithilfe des Parameters $q \in \{0,1,2,3,4,5\}$ kann die Anzahl $2^q$ der Auswertungsstellen für die Gauss-Quadratur gewählt werden. Die beiden Parameter $\zeta_0,\zeta_1 \in \R^+$ werden in den Zulässigkeitsbedingungen verwendet. Der Parameter \lstinline!typ! $\in \{1,2,3\}$ steht für die Art der Berechnung.
+aufgerufen werden kann. Hierbei steht wie in der Datenstruktur beschrieben \lstinline!coo! für die Koordinatenmatrix und \lstinline!ele! für die Elementmatrix. Weiterhin sei \lstinline!zeta! ein Array mit den Einträgen \lstinline!zeta! $=(q,\zeta_Q,\zeta_E)$. Mithilfe des Parameters $q \in \{0,1,2,3,4,5\}$ kann die Anzahl $2^q$ der Auswertungsstellen für die Gauss-Quadratur gewählt werden. Die beiden Parameter $\zeta_Q,\zeta_E \in \R^+$ werden in den Zulässigkeitsbedingungen verwendet. Der Parameter \lstinline!typ! $\in \{1,2,3\}$ steht für die Art der Berechnung.
\noindent
Insgesamt wurden vier Arten der Berechnung implementiert.
\end{itemize}
\item In der vollen Quadratur für zulässige Elemente werden
\begin{itemize}
- \item $\zeta_0$-zulässige Elemente mit Gauss-Quadratur über alle Integrale
+ \item $\zeta_Q$-zulässige Elemente mit Gauss-Quadratur über alle Integrale
\item und alle anderen analytisch mittels Stammfunktion berechnet.
\end{itemize}
\item Bei der geeigneten Quadratur für zulässige Elemente werden
\begin{itemize}
- \item $\zeta_1$-zulässige Elemente mit Gauss-Quadratur über ein Element,
- \item $\zeta_0$-zulässige Elemente mit Gauss-Quadratur über alle Integrale
+ \item $\zeta_E$-zulässige Elemente mit Gauss-Quadratur über ein Element,
+ \item $\zeta_Q$-zulässige Elemente mit Gauss-Quadratur über alle Integrale
\item und alle anderen analytisch mittels Stammfunktion berechnet.
\end{itemize}
\item In der Quadratur über ein Element für zulässige Elemente werden
\begin{itemize}
- \item $\zeta_0$-zulässige Elemente mit Gauss-Quadratur über ein Element
+ \item $\zeta_E$-zulässige Elemente mit Gauss-Quadratur über ein Element
\item und alle anderen analytisch mittels Stammfunktion berechnet.
\end{itemize}
\end{enumerate}
}
\subsection{Fehlerschätzer}
-Mithilfe des Galerkin-Verfahrens berechnen wir nur eine approximative Lösung, für die wir eine Aussage über die Genauigkeit der Lösung treffen wollen. Da wir für den Fehler $\enorm{\phi - \phi_{\ell}}$ nur die Galerkin-Lösung $\phi_{\ell}$ kennen und die exakte Lösung $\phi$ unbekannt ist, werden wir im Folgenden verschiedene Fehlerschätzer betrachten.
+Mithilfe des Galerkin-Verfahrens berechnen wir nur eine approximative Lösung von \eqref{math:slp:gls}, für die wir eine Aussage über die Genauigkeit der Lösung treffen wollen. Da wir für den Fehler $\enorm{\phi - \phi_{\ell}}$ nur die Galerkin-Lösung $\phi_{\ell}$ kennen und die exakte Lösung $\phi$ unbekannt ist, werden wir im Folgenden verschiedene Fehlerschätzer betrachten.
% In diesem Abschnitt definieren wir die a-posteriori Fehlerschätzer, die wir zur Steuerung des adaptiven Algorithmus einsetzen werden.
Im Wesentlichen werden wir hierzu die $h-h/2$ Strategie aus \cite{fer:errbem} verwenden.
-Im Folgenden bezeichnen wir mit $\hat \T_{\ell}$ das Gitter welches entsteht, wenn das Gitter $\T_{\ell}$ uniform, also entlang aller Kanten geteilt wird. Weiterhin bezeichne $\phi$ die exakte Lösung des Galerkin-Verfahrens und $\phi_{\ell}$ die Lösung zum Gitter $\T_{\ell}$, sowie $\hat \phi_{\ell}$ die Lösung zum uniformen Gitter $\hat \T_{\ell}$.
-\begin{defi}Es bezeichne $\phi$ die exakte Lösung von Formel \todo{ref}, $\phi_{\ell}$ die Galerkin-Lösung auf dem Gitter $\T_{\ell}$ und $\hat \phi_{\ell}$ die Lösung auf dem uniform verfeinerten Gitter $\hat \T_{\ell}$. Dann gilt, der Schätzer
+Im Folgenden bezeichnen wir mit $\hat \T_{\ell}$ das Gitter welches entsteht, wenn das Gitter $\T_{\ell}$ uniform, also entlang aller Kanten geteilt wird. Weiterhin bezeichne $\phi$ die exakte Lösung des Galerkin-Verfahrens und $\phi_{\ell}$ die Galerkin-Lösung zum Gitter $\T_{\ell}$, sowie $\hat \phi_{\ell}$ die Galerkin-Lösung zum uniformen Gitter $\hat \T_{\ell}$.
+\begin{defi}Es bezeichne $\phi$ die exakte Lösung von \eqref{math:slp:gls}, $\phi_{\ell}$ die Galerkin-Lösung von \eqref{math:slp:gls:galerkin} auf dem Gitter $\T_{\ell}$ und $\hat \phi_{\ell}$ die Galerkin-Lösung auf dem uniform verfeinerten Gitter $\hat \T_{\ell}$. Dann gilt, der Schätzer
\begin{align*}
\eta_{\ell} &:= \enorm{\hat \phi_{\ell} - \phi_{\ell}}
\end{align*}
\end{align*}
und unter der Saturationsannahme
\begin{align*}
-\norm{\phi -\hat \phi_{\ell}} &\leq C_{sat} \cdot \norm{\phi - \phi_{\ell}} \quad \text{mit }0 < C_{sat} < 1
+\norm{\phi -\hat \phi_{\ell}} &\leq C_{sat} \norm{\phi - \phi_{\ell}} \quad \text{mit }0 < C_{sat} < 1
\end{align*}
zuverlässig
\begin{align*}
\subsection{Markieren}
Im Adaptiven Algorithmus werden wir die Elemente abhängig vom Fehlerschätzer $\tilde \mu$ verfeinern. Dazu wählen wir mithilfe der Dörfler-Markierung \cite{dor:adapt} eine Teilmenge aus.
-\begin{defi}[Dörfler-Markierung]
+\begin{defi}[Dörfler-Markierung]\label{thm:def:mark}
Bestimme für feste Konstante $\theta \in (0,1]$ die Menge $M_{\ell} \subseteq T_{\ell}$ mit minimaler Kardinalität
\begin{align*}
\theta \sum_{T\in T_{\ell}} \tilde \mu_{\ell}(T)^2 &\leq \sum_{ T\in M_{\ell}} \tilde \mu_{\ell}(T)^2.
Weiterhin sei $marked_j = 1$ für alle $T_j \in \T_{\ell} \backslash \tilde M_{\ell}$.
\end{defi}
Die Funktion
-\begin{align*}
-marked = mark(xF2S, tmu, theta, nu); \quad \text{mit } xF2S := x_{fine}[F2S]
-\end{align*}
-implementiert die Definitionen zum Bestimmen der Markierung.
+\begin{lstlisting}[language=M,numbers=none]
+marked = mark(xF2S, tmu, theta, nu);
+\end{lstlisting}
+implementiert die Definitionen zum Bestimmen der Markierung. Hierbei seien die Einträge der Matrix $xF2S \in \R^{N\times4}$ geben durch $(xF2S)_{ij} = \langle\hat\phi_{\ell},\chi_{(F2S_{ij})}\rangle$.
\begin{enumerate}
\renewcommand{\theenumi}{(\roman{enumi})}
\item Verfeinere $T_{\ell}$ um $\hat T_{\ell}$ zu erhalten \label{alg:adapt:begin}
-\item Berechne die Galerkinlösung $\hat \phi_{\ell} \in P^0(\hat T_{\ell})$
+\item Berechne die Galerkin-Lösung $\hat \phi_{\ell} \in P^0(\hat T_{\ell})$ von \eqref{math:slp:gls:galerkin}
\item Berechne Fehlerschätzer $\tilde \mu_{i} := \norm{\varrho^{\ell/2}(\hat \phi_{\ell} - \Pi_{\ell} \hat \phi_{\ell} )}$
-\item Wähle $M_{\ell} \subseteq T_{\ell}$ mit minimaler Kardinalität, so dass
-\begin{align*}
-\theta \sum_{T\in \T_{\ell}} \tilde\mu_{\ell}^2 & \leq \sum_{T\in M_{\ell}} \tilde\mu_{\ell}^2
-\end{align*}
-\item Verfeinere mindestens die Markierten Elemente $M_{\ell}$ um $\T_{\ell+1}$ zu erhalten
+\item Wähle Teilmenge $M_{\ell} \subseteq T_{\ell}$ wie in Definition \ref{thm:def:mark}
+\item Verfeinere mindestens die Markierten Elemente $M_{\ell}$ durch Algorithmus \ref{alg:refine} um $\T_{\ell+1}$ zu erhalten
\item $\ell \mapsto \ell+1$, gehe zu \ref{alg:adapt:begin}
\end{enumerate}
\end{alg}
\begin{align*}
\{ (0,0,0), (1,0,0), (1,1,0), (0,1,0) \}.
\end{align*}
-Die Lösung $\phi$ kann nicht exakt bestimmt werden, weist aber an den vier Eckpunkten Singularitäten gegen $\infty$ auf. Wir werden aber trotzdem die Energienorm als Referenzlösung mithilfe des {\it Aitken'schen} $\varDelta^2$-Verfahren bestimmen. Das Verfahren dient zum Beschleunigen der Konvergenz von Folgen. Damit erhalten wir für die exakte Lösung $\phi$ die Energienorm $\enorm{\phi}^2 = 4.609193$.
+Die Lösung $\phi$ kann nicht exakt bestimmt werden, weist aber an den vier Eckpunkten Singularitäten auf. Wir werden aber trotzdem die Energienorm als Referenzlösung mithilfe des {\it Aitken'schen} $\varDelta^2$-Verfahren bestimmen. Das Verfahren dient zum Beschleunigen der Konvergenz von Folgen. Damit schätzen wir für die exakte Lösung $\phi$ die Energienorm $\enorm{\phi}^2 \approx 4.609193$.
\begin{figure}[ht]
\centering
\subfloat[Startnetz \label{fig:mesh:2DQuad:start}]{\includegraphics[width=0.5\textwidth]{fig/exmpl_2DQuad_ref}}
Um auch die Stabilität der drei Strategien untersuchen zu können, sehen wir in der Abbildung \ref{fig:2DQuad:verfeinern:cond} die Konditionszahlen der $V_{\ell}$ Matrix in Abhängigkeit der Elementanzahl.
\noindent
-Weiterhin können wir in Abbildung \ref{fig:2DQuad:verfeinern:time} die benötigte Zeit pro Berechnungsschritt ablesen. Hierbei fällt auf, dass die Wahl der Strategie keinen Einfluss auf die benötigte Zeit hat, sondern nur die Anzahl der Elemente. Für die Berechnung mit 3000 Elementen benötigen alle drei Strategien etwa $10^4s \approx 2h45m$.
+Weiterhin können wir in Abbildung \ref{fig:2DQuad:verfeinern:time} die benötigte Zeit pro Berechnungsschritt ablesen. Hierbei fällt auf, dass die Wahl der Strategie keinen Einfluss auf die benötigte Zeit hat, sondern nur die Anzahl der Elemente. Für die Berechnung mit 3000 Elementen benötigen alle drei Strategien etwa $10^4$ Sekunden, was fast drei Stunden entspricht.
\noindent
-Diese Ergebnisse Zeigen also, dass die "`anisotrope"' Strategie die beste Konvergenzrate aufweist, wir dafür jedoch eine schlechtere Konditionsrate in kauf nehmen müssen. Welches letztendlich auf die größen und formen der Elemente zurückzuführen ist. An dieser Stelle wollen wir noch zusätzlich Abbildung \ref{fig:mesh:2DQuad:steps} betrachten, welche das "`anisotrop"' verfeinerte Netz nach 12 Schritten darstellt. Denn hier erkennen wir sehr gut, dass diese Strategie das Netz insbesondere an den Singularitäten verfeinert.
+Diese Ergebnisse Zeigen also, dass die "`anisotrope"' Strategie die beste Konvergenzrate aufweist, wir dafür jedoch eine schlechtere Konditionszahl der Matrix in kauf nehmen müssen. Dies ist letztendlich auf die Größe und Forme der Elemente zurückzuführen. An dieser Stelle wollen wir noch zusätzlich Abbildung \ref{fig:mesh:2DQuad:steps} betrachten, welche das "`anisotrop"' verfeinerte Netz nach 12 Schritten darstellt. Denn hier erkennen wir sehr gut, dass diese Strategie das Netz insbesondere an den Singularitäten verfeinert.
\noindent
Ziel wird es nun sein die Instabilitäten, die bei der analytischen Berechnung auftreten durch Quadratur zu vermeiden. Dafür werden wir vorher noch Berechnungen mit verschiedenen Quadraturgraden genauer untersuchen.
\subsubsection{Vergleich verschiedener Quadraturgrade}
-Bei der folgenden Berechnung werden wir alle Integrale von $\zeta_Q$-zulässigen Elemente durch die vorgestellte Gauss-Quadratur approximieren. Hierbei werden wir jeweils 1, 2, 4 oder 8 Auswertungsstellen verwenden. Alle vier Berechnungsarten werden auf dem selben Netz ausgeführt um die Ergebnisse nicht durch die Wahl des Netzes zu beeinflusse. Zum Berechnen des $\tilde \mu_{\ell}$-Schätzers, welcher die Verfeinerung steuert, werden wir in jedem Schritt die Lösung der Quadratur mit 8 Auswertungsstellen verwenden.
+Bei der folgenden Berechnung werden wir alle Integrale von $\zeta_Q$-zulässigen Elemente durch die vorgestellte Gauss-Quadratur approximieren. Hierbei werden wir jeweils 1, 2, 4 oder 8 Auswertungsstellen verwenden. Alle vier Berechnungsarten werden auf den selben Netzen ausgeführt um die Ergebnisse nicht durch die Wahl des Netzes zu beeinflusse. Zum Berechnen des $\tilde \mu_{\ell}$-Schätzers, welcher die Verfeinerung steuert, werden wir in jedem Schritt die Lösung der Quadratur mit 8 Auswertungsstellen verwenden.
\begin{figure}[ht]
In Abbildung \ref{fig:2DQuad:quad:err} haben wir die Ergebnisse der Fehler und Fehlerschätzer für die verschiedenen Quadraturgrade dargestellt. Wir beobachten, dass eine sowie zwei Auswertungsstellen, dargestellt durch die Linien in \figLineA[] und \figLineB[] instabil werden. Für die Quadraturgrade vier und acht dargestellt durch die Linien in \figLineC[] und \figLineD[], erkennen wir, dass die Berechnungen stabil bleiben und das auch für Elementanzahlen bei denen die analytische Berechnung (vergleiche Abbildung \ref{fig:2DQuad:verfeinern}) versagt.
\noindent
-Da wir möglichst rechen ökonomisch arbeiten wollen, haben wir auch die Berechnungszeiten in Abbildung \ref{fig:2DQuad:quad:time} untersucht. Wie sich leicht erkennen lässt sind die Berechnungszeiten für die verschiedenen Quadraturgrade etwa äquivalent. Lediglich die Quadratur mit 8 Auswertungsstellen benötigt für kleine Elementanzahlen etwas länger, wobei sie für große Anzahlen jedoch fast gleich ist. Daraus lässt sich ablesen, dass die Wahl des Quadraturgrades für große Netze kaum einen Einfluss auf die Berechnungszeit nimmt.
+Da wir möglichst ökonomisch arbeiten wollen, haben wir auch die Berechnungszeiten in Abbildung \ref{fig:2DQuad:quad:time} untersucht. Wie sich leicht erkennen lässt sind die Berechnungszeiten für die verschiedenen Quadraturgrade etwa äquivalent. Lediglich die Quadratur mit 8 Auswertungsstellen benötigt für kleine Elementanzahlen etwas länger, wobei sie für große Anzahlen jedoch fast gleich ist. Daraus lässt sich ablesen, dass die Wahl des Quadraturgrades für große Netze kaum einen Einfluss auf die Berechnungszeit nimmt.
\noindent
Aufgrund dieser Ergebnisse werden wir für die Folgenden Berechnungen einen Quadraturgrad von 4 wählen. Denn wir wollen zum einen die Berechnungszeiten auch für kleine Netze gering halten und zum anderen die Stabilität der Berechnungen sicher stellen.
% \begin{align*}
% \{ (0,0,0), (1,0,0), (1,1,0), (0,1,0) \}.
% \end{align*}
-Die Lösung $\phi$ strebt in diesem Fall an allen außen liegenden Eckpunkten gegen $+\infty$ und an dem Innen liegenden Eckpunkt $(0, 0, 0)$ gegen $-\infty$.
-Auch hier wurde die Energienorm $\enorm{\phi}^2 = 16.2265$ der exakten Lösung $\phi$ mithilfe des {\it Aitken'schen} $\varDelta^2$-Verfahrens bestimmt.
+Die Lösung $\phi$ strebt in diesem Fall an allen außen liegenden Eckpunkten gegen $+\infty$ und an dem innen liegenden Eckpunkt $(0, 0, 0)$ gegen $-\infty$.
+Auch hier wurde die Energienorm $\enorm{\phi}^2 \approx 16.2265$ der exakten Lösung $\phi$ mithilfe des {\it Aitken'schen} $\varDelta^2$-Verfahrens geschätzt.
\begin{figure}[ht]
\centering
\subfloat[Startnetz \label{fig:mesh:3DFichCube:start}]{\includegraphics[width=0.5\textwidth]{fig/exmpl_3DFichCube_ref}}
\item plot
\end{enumerate}
}
- \subsection{C++}
+ \subsection[C++]{\Cpp}
+ An dieser Stelle werden wir den Quellcode für die Berechnung der Matrix $A$ aus Kapitel \ref{sec:semi} und \ref{sec:analyt} vorstellen.
\subsubsection{mex\_build\_V.cpp}\label{code:mex_build_V}
- \lstinputlisting[language=C++]{../src/mex_build_V.cpp}
+ Diese Datei enthält die Implementierung für die Berechnung der Matrix $A$ in {\sc MEX}. Sie benötigt zum Kompilieren die Dateien \lstinline!slpRectangle.hpp!, \lstinline!slpRectangle.cpp! und \lstinline!gauss.hpp!.\\
+ Innerhalb von \Matlab~kann die Funktion aufgerufen werden mit
+\begin{lstlisting}[language=M,numbers=none]
+[A] = mex_build_V(coo,ele,zeta,typ);
+\end{lstlisting}
+ Hierbei ist \lstinline!coo! $\in \R^{M\times3}$ die Koordinatenmatrix aus Kapitel \ref{sec:implement:daten} und \lstinline!ele! $\in \N^{N\times4}$ die zugehörige Elementmatrix in der die Zeilenindizes der Koordinaten \lstinline!coo! stehen. Eine Zeile der Koordinatenmatrix entspricht einer Koordinate und eine Zeile aus der Elementmatrix steht für die Indizes der Koordinaten eines Elements. Weiterhin steht \lstinline!zeta! $ = (p, \zeta_Q, \zeta_E)$ für den Koordinatenvektor. Mithilfe des Parameters $q \in \{0,1,2,3,4,5\}$ kann die Anzahl $2^q$ der Auswertungsstellen für die Gauss-Quadratur gewählt werden. Die beiden Parameter $\zeta_Q,\zeta_E \in \R^+$ werden in den Zulässigkeitsbedingungen verwendet. Der Parameter \lstinline!typ! $\in \{1,2,3\}$ steht für die Art der Berechnung aus Kapitel \ref{sec:implement:A}. Zurück gegeben wird die zur Parametrisierung gehörende Matrix $A \in\R^{N\times N}$.
+ \lstinputlisting[language=C++,firstline=29,firstnumber=29]{../src/mex_build_V.cpp}
\subsubsection{slpRectangle.hpp}\label{code:slpRectangle.hpp}
- \lstinputlisting[language=C++]{../src/slpRectangle.hpp}
+ Diese Datei stellt die Funktionen zur Berechnung eines Matrixeintrages $A_{jk}$ bereit, welcher dann für parallele Elemente über
+\begin{lstlisting}[language=C++,numbers=none]
+double cParO1(b,d,t,v,d1,d2,d3,zeta*);
+\end{lstlisting}
+ berechnet werden kann. \lstinline!b,d! entspricht den Seitenlängen des Elements $T_j$ und \lstinline!t,v! denen des Elements $T_k$ und $\bs \delta =$\lstinline!(d1,d2,d3)! dem Abstandsvektor, wie in Kapitel \ref{sec:analyt:int}. Hierbei ist zu beachten, dass \lstinline!b,t! in der selben Dimension liegen wie auch \lstinline!d1! und wie auch \lstinline!d,v,d2! in der selben Dimension liegen. Die Zulässigkeitsbedingungen werden durch \lstinline!zeta*! = $(\zeta_Q,\zeta_E)$ bestimmt. Für Ortogonale Elemente steht analog die Funktion
+\begin{lstlisting}[language=C++,numbers=none]
+double cOrtO1(b,d,t,v,d1,d2,d3,zeta*);
+\end{lstlisting}
+zur Verfügung. Hier liegen jedoch, wie in Kapitel \ref{sec:analyt:int} für orthogonale Elemente vorgestellt, \lstinline!b,d1! und \lstinline!d,t,d2! und \lstinline!v,d3! in jeweils der selben Dimension.
+ \lstinputlisting[language=C++,firstline=13,firstnumber=13]{../src/slpRectangle.hpp}
\subsubsection{slpRectangle.cpp}\label{code:slpRectangle.cpp}
- \lstinputlisting[language=C++]{../src/slpRectangle.cpp}
+ Diese Datei implementiert die Funktionen aus \lstinline!slpRectangle.hpp!. Hierzu wurden die vorgestellten Funktionen aus Kapitel \ref{sec:analyt} umgesetzt, wie auch die zugehörigen Quadraturen.
+ \lstinputlisting[language=C++,firstline=13,firstnumber=13]{../src/slpRectangle.cpp}
\subsubsection{gauss.hpp}\label{code:gauss}
- \lstinputlisting[language=C++]{../src/gauss.hpp}
- \subsection{Matlab}
+ Diese Datei wurde automatisch von der \Matlab-Funktion \lstinline!export_gauss.m! erstellt.
+ Durch den Befehl
+\begin{lstlisting}[language=M,numbers=none]
+export_gauss(start,stop,steps,file);
+\end{lstlisting}
+ wird die Datei \lstinline!file! angelegt, welche dann in \Cpp verwendet werden kann. In ihr sind die Gauss-Quadraturpunkte und Gewichte für das Intervall $[start, stop]$ gespeichert. \lstinline!steps! steht hierbei für einen Vektor, welcher die zur Verfügung stehenden Grade steuert.
+
+ \noindent
+ In \Cpp~kann die Quadratur dann wie im folgenden Beispiel mit Quadraturgrad $2^3$ benutzt werden.
+\begin{lstlisting}[language=C++]
+double sol;
+for(int i=0;i<GAUSS_SIZE[3];++i)
+ sol += sin(GAUSS_NODES[3][i].n) * GAUSS_NODES[3][i].w;
+\end{lstlisting}
+
+ \noindent
+ Da die Datei fast nur Quadraturpunkte und Gewichte enthält wollen wir hier nur den Anfang vorstellen.
+ \lstinputlisting[language=C++,firstline=26,firstnumber=26,lastline=41]{../src/gauss.hpp}
+
+ \subsection[Matlab]{\Matlab}
+ Hier wollen wir alle \Matlab-Funktionen für den adaptiven Algorithmus \ref{alg:adapt} vorstellen. Funktionen die lediglich zum Erstellen von Grafiken und zum Testen dienen werden wir hier nicht zeigen.
\subsubsection{compute.m}\label{code:compute}
- \lstinputlisting[language=M]{../src/compute.m}
+ Diese Funktion Implementiert im Wesentlichen den Algorithmus \ref{alg:adapt}. Nachdem das Startnetz aus der Datei \lstinline!file! geladen wurde, werden \lstinline!times! $\in\N$ Verfeinerungsschritte durchgeführt. Sollte \lstinline!times! größer als 40 sein, wird nicht nach 40 Schritten abgebrochen sondern nachdem die Elementanzahl \lstinline!times! vom Netz $\T_{\ell}$ erreicht wurde. Die Berechnungsarten werden über den Vektor \lstinline!typ! $\in\{1,2,3,4\}^n$ bestimmt. Hierbei werden in jedem Verfeinerungsschritt alle $n$ Berechnungsarten aus \lstinline!typ! auf dem selben Netz durchgeführt. \lstinline!zeta! ist ein Array welcher zu den $n$ Berechnungsarten den entsprechenden $\zeta=(p,\zeta_Q,\zeta_E)$ Vektor enthält. Zum steuern des Verfeinernalgorithmus \ref{alg:refine} dienen die Parameter \lstinline!theta,nu! $\in[0,1]$ und über den Parameter \lstinline!vcon! $\in\{0,1\}$ kann die Vorkonditionierung der $A$ Matrix eingeschaltet werden. Das neue Netz $\T_{\ell+1}$ wird durch die letzte Berechnungsart bestimmt. Zurück gegeben wird die Matrix \lstinline!data!, in der alle wichtigen Auswertungen gespeichert sind.
+ Die folgende Zeile
+\begin{lstlisting}[language=M,numbers=none]
+compute('exmpl_2DQuad',500,{[0 0 0] [3 2 2]},[1 2],0.5,0.5,0);
+\end{lstlisting}
+führt eine adaptiv anisotrope ($\theta=0.5,\nu=0.5$) voll Analytische und volle Quadratur, für $\zeta_Q=2$-Zulässige Elemente, (\lstinline!typ=[1 2]!) Berechnung bis \lstinline!times=500! Elemente auf einem Quadrat durch. Hierbei wird eine Gauss-Quadratur vom Grad 8 verwendet, keine Vorkonditionierung und zum Markieren der zu verfeinernden Elemente die Galerkin-Lösung zum Typ "`volle Quadratur"'. Gelöst wird das Problem aus \eqref{math:bsp:Quad:gls}.
+\begin{lstlisting}[language=M]
+function [data er fileo] = compute(file,times,zeta,typ,theta,nu,vcon)
+\end{lstlisting}
+ \lstinputlisting[language=M,firstline=19,firstnumber=19]{../src/compute.m}
\subsubsection{refineQuad.m}\label{code:refineQuad}
- \lstinputlisting[language=M]{../src/refineQuad.m}
- \subsubsection{areaQuad.m}\label{code:areaQuad}
- \lstinputlisting[language=M]{../src/areaQuad.m}
+ Diese Funktion implementiert den Verfeinern Algorithmus \ref{alg:refine}. Übergeben werden vom Netz $\T_{\ell}$, die Koordinatenmatrix \lstinline!coordinates! $\in\R^{M\times3}$, die Elementmatrix \lstinline!elements! $\in\{1,\ldots,M\}^{N\times4}$, die Nachbarschaftsrelationen \lstinline!neigh! $\in\{0,1,\ldots,N\}^{N\times8}$. Weiterhin wird auch die Seitenmatrix \lstinline!sites! $\in\N^{N\times2}$ übergeben. Die $j$te Zeile $(a_{j},b_{j})$ entspricht der Anzahl der Halbierungen der Seitenlängen $a,b$ vom Element $T_j$, gelten muss hierbei $a = 2^{a_j}$ und $b = 2^{b_j}$. Außerdem entspricht \lstinline!typ! $\in\{1,\ldots,5\}^N$ der Markierung. Bei der Rückgabe des verfeinerten Netzes $\tilde{\T_{\ell}}$ werden zusätzlich die VaterSohn Beziehungen \lstinline!f2s! zurückgegeben.
+\begin{lstlisting}[language=M]
+function [coo,ele,nei,f2s,sit,err] = refineQuad(coordinates,elements,neigh,sites,typ)
+\end{lstlisting}
+ \lstinputlisting[language=M,firstline=18,firstnumber=18]{../src/refineQuad.m}
+% \subsubsection{areaQuad.m}\label{code:areaQuad}
+% Diese Funktion berechnet den Flächeninhalt \lstinline!area! $\in\R^N$ und die exakten Seitenlängen \lstinline!site! $\in\R^{N\times2}$ aller Elemente anhand der Matrix \lstinline!site! $\in \N^{N\times2}$.
+% \begin{lstlisting}[language=M]
+% function [area, sites] = areaQuad(sites)
+% \end{lstlisting}
+% \lstinputlisting[language=M,firstline=30,firstnumber=30]{../src/areaQuad.m}
\subsubsection{mark.m}\label{code:mark}
- \lstinputlisting[language=M]{../src/mark.m}
+ Diese Funktion Implementiert die Definition \ref{thm:def:mark}.
+\begin{lstlisting}[language=M]
+function REF = mark(xF2S,ind,theta,nu)
+\end{lstlisting}
+ \lstinputlisting[language=M,firstline=9,firstnumber=9]{../src/mark.m}
\clearpage
\bibliographystyle{gerabbrv}