]> git.leopard-lacewing.eu Git - bacc.git/commitdiff
[doc] viel neuer Text
authorPeter Schaefer <schaeferpm@gmail.com>
Sun, 21 Apr 2013 14:13:01 +0000 (16:13 +0200)
committerPeter Schaefer <schaeferpm@gmail.com>
Sun, 21 Apr 2013 14:13:01 +0000 (16:13 +0200)
[src] count versucht in mex_build einzuführen
[src] areaQuad wieder auf coo,ele umgestellt und in compute direkt berechnet
[src] Kommentare mehr ausgemistet

doc/doc.pdf
doc/doc.tex
src/areaQuad.m
src/compute.m
src/export_gauss.m
src/mex_build_V.cpp
src/refineQuad.m
src/slpRectangle.cpp
src/slpRectangle.hpp
src/test_sol.m

index a9dcc8eb45b9c6da8d7e0b24be087b4741747b23..1912e050b74725a26c0ffee70a0535e24cbbb82d 100644 (file)
Binary files a/doc/doc.pdf and b/doc/doc.pdf differ
index f313739194d1c0aa64aa448139ce264abef65c65..99d9e101ce8977010f7bcb4366b8ae6b106cf0ed 100644 (file)
@@ -315,26 +315,26 @@ Sei nun $\phi$ die eindeutige Lösung von \eqref{math:slp:gls} und bezeichne $g$
 
 \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}
@@ -611,7 +611,7 @@ wobei die Lagrange-Polynome $L_j$ definiert sind durch
 \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*}
@@ -1553,7 +1553,7 @@ In diesem Abschnitt wollen wir uns mit der analytischen Berechnung des Integrals
   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
@@ -1579,7 +1579,7 @@ g(-3/2;y;x;\lambda) &= \frac{y-x}{\lambda^2}\{(y-x)^2+\lambda^2\}^{-1/2}
 \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\\
@@ -1606,8 +1606,8 @@ Für alle weiteren relevanten Fälle $p = -3/2 +k$ mit $k \in \N$ können wir fo
 \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
@@ -1651,7 +1651,7 @@ Analog können wir orthogonal liegende Elemente mit Definition \ref{thm:def:T} s
 \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}\\
@@ -1685,7 +1685,7 @@ Dann gilt nach \cite{mai:3dbem}
 
 
 \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} 
@@ -1793,7 +1793,7 @@ Hierzu werden wir kurz Festhalten wie die Diskretisierung des Netzes aus Kapitel
 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}
@@ -1803,16 +1803,17 @@ Die Elemente $\T_{\ell}$ werden wir ebenfalls Zeilenweise in einer $N \times 4$
 \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
 
@@ -1833,8 +1834,10 @@ Offensichtlich ist $i \notin N_i$. Wir wollen uns aber noch genauer eine geeigne
 
 
 \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}
 
@@ -1844,10 +1847,10 @@ Da wir im weiteren Verlauf sowohl adaptive also auch anisotrope Netzverfeinerung
 \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.
@@ -1856,7 +1859,7 @@ Damit jedem Element $T_i$ eine Teilungsart zugeordnet werden kann, haben wir ein
 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);
@@ -1873,7 +1876,7 @@ Da wir später den Fehlerschätzer berechnen wollen, ist es wichtig sich zu jede
 \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}.
@@ -1899,14 +1902,14 @@ Die Implementierung des Integrals \eqref{math:imple:Ajk} wurde in \lstinline!mex
 \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.
@@ -1917,18 +1920,18 @@ 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}
@@ -1960,11 +1963,11 @@ Die Implementierung der Berechnung eines Eintrags $A_{jk}$ befindet sich in der
 }
 
 \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*}
@@ -1974,7 +1977,7 @@ ist effizient
 \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*}
@@ -2033,7 +2036,7 @@ Dadurch gilt auf isotropen Netzen:
 
 \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.
@@ -2075,10 +2078,10 @@ marked_j :=
 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$.
 
 
 
@@ -2101,13 +2104,10 @@ Mithilfe der oben Definierten Funktionen ist es uns nun möglich den Ablauf der
 \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}
@@ -2124,7 +2124,7 @@ mit dem Startnetz aus Abbildung \ref{fig:mesh:2DQuad:start}, einem Quadrat in de
 \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}}
@@ -2204,16 +2204,16 @@ Anhand der \figLineC[] farbenen Linien, also der "`anisotropen"' Strategie beoba
 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]
 
@@ -2232,7 +2232,7 @@ Bei der folgenden Berechnung werden wir alle Integrale von $\zeta_Q$-zulässigen
 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.
@@ -2266,8 +2266,8 @@ Der Fischer Würfel ist gegeben durch einen $[-1, 1]^3$ Würfel aus dem ein $[-1
 % \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}}
@@ -2372,24 +2372,81 @@ Die wichtigsten Funktionen
 \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}
index 1f9652dde453e3a62b59428c18f9ac9ad9cf629f..c95ba2fe1a77ce2a92e26e8616a7eedcc88bc10a 100644 (file)
@@ -1,4 +1,4 @@
-function [area, sites] = areaQuad(sites)
+function [area, vector, sites] = areaQuad(coordinates, elements,varargin)
 %
 % [area, vector, sites] = area(coordinates, elements)
 %
@@ -8,27 +8,27 @@ function [area, sites] = areaQuad(sites)
 %
 % P.Schaefer
 
-    vector = zeros(size(elements,1),3);
-    area = zeros(size(elements,1),1);
-    sites = zeros(size(elements,1),2);
+    vector = zeros(size(elements,1),3);
+    area = zeros(size(elements,1),1);
+    sites = zeros(size(elements,1),2);
     
 %%  Flaecheninhalt und Normalen berechnen
 
-%     for i = 1:size(elements,1)
-%         % normalized Vector on every triangle
-%         tri = elements(i,:);
-%         a = (coordinates(tri(2),:)-coordinates(tri(1),:));
-%         b = (coordinates(tri(4),:)-coordinates(tri(1),:));
-%         N = cross(a',b');
-%         
-% %         area(i) = norm(N);
-%         
-%         N = N/norm(N);
-%         vector(i,:) = N;
-        area =  2.^-sum(sites,2);
-        sites = 2.^-sites ;
+    for i = 1:size(elements,1)
+        % normalized Vector on every triangle
+        tri = elements(i,:);
+        a = (coordinates(tri(2),:)-coordinates(tri(1),:));
+        b = (coordinates(tri(4),:)-coordinates(tri(1),:));
+        N = cross(a',b');
         
-%     end
+%         area(i) = norm(N);
+        
+        N = N/norm(N);
+        vector(i,:) = N;
+        area(i) = norm(a) * norm(b);
+        sites(i,:) = [norm(a) norm(b)];
+        
+    end
 
 end
 
index 8e67a8050d390705bbae7c88c2d685e9c578bc86..068eb30bcff88594b246d1ea6267b7e3a712c6aa 100644 (file)
@@ -46,9 +46,9 @@ for j = 1:times
     =refineQuad(coordinates,elements,neigh,sites,2);
 
   %Flaecheninhalte Berechnen (rhs)
-  b_fine = areaQuad(sit_fine);
+  b_fine = 2.^-sum(sit_fine,2);
   
-  b = areaQuad(sites);
+  b = 2.^-sum(sites,2);
   hmin = 2.^-max(sites,[],2);
   hmax = 2.^-min(sites,[],2);
 
@@ -211,13 +211,13 @@ for j = 1:times
     = refineQuad(coordinates,elements,neigh,sites,marked);
   
   %Vater Sohn test
-  assert(sum(areaQuad(old_S))==sum(areaQuad(sites)),...
+  assert(sum(2.^-sum(old_S,2))==sum(2.^-sum(sites,2)),...
     'Gesamtinhalt Fehlerhaft')
   
   e = 0;
   for k = 1:size(data,1)
-    e = e +  abs(areaQuad(old_S(k,:))...
-      - sum(areaQuad(sites(unique(f(k,:)'),:))));
+    e = e +  abs(2.^-sum(old_S(k,:))...
+      - sum(2.^-sum(sites(unique(f(k,:)'),:),2)));
   end
   assert(e==0, 'Einzelne Inhalte Fehlerhaft')
   
@@ -250,20 +250,6 @@ for j = 1:times
   save (['meshSave/' fileo int2str(size(data,1))]...
     , 'coordinates', 'elements','neigh','data', 'sites')
   
-  %Alle Relevanten zwischenInformationen Speichern
-%   out = '_';
-%   if(length(varargin)~=0)
-%      out = varargin{1};
-%    end
-%   typeN = int2str(type);
-%  save (['meshSave/fine_' out typeN(typeN~=' ')...
-%      '_' int2str(size(data,1))],...
-%      'coo_fine', 'ele_fine','neigh_fine','f2s','data',...
-%      'save_A','save_x','save_A_fine','save_x_fine','b','b_fine')
-
-%   clear 'coo_fine' 'ele_fine' 'neigh_fine' 'f2s' 
-%   plotShape(G_C,G_E,'');
-
 end
 end
 
index 681e4522593fbad744686e22ab1d99e228b93e59..91fb9d968780013d48f9efe39be6329815cc8d7c 100644 (file)
@@ -60,8 +60,8 @@ char(10) ...
 char(10)];
 
 preci = 21;
-str2 =['gauss const GAUSS_NODES[][' num2str(max(steps)) '] = {' char(10)];
-str1 =['double const GAUSS_SIZE[] = { '];
+str2 =['gauss GAUSS_NODES[][' num2str(max(steps)) '] = {' char(10)];
+str1 =['double GAUSS_SIZE[] = { '];
 
 for i = steps(1:end)
   str1 = [str1 num2str(i) ', '];
index 8ae813b37d808e977822a50d166077c03b2a60c5..4a8867c7599ce168ba933f37eb9431fa78938ecb 100644 (file)
@@ -29,7 +29,7 @@
 /*\r
  * wenn gesetzt wird paralleles Rechnen benutzt\r
  */\r
-//#define PARALLEL\r
+#define PARALLEL\r
 #include <cmath>\r
 #include <mex.h>\r
 \r
@@ -214,6 +214,7 @@ double inline distT_ort(double b, double d, double t, double v, double d1,
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
 \r
        int i, j, k; //Schleifenindizes\r
+       bool count = false;\r
 \r
        /*\r
         * Sicherheitsueberpruefung der Parameter\r
@@ -221,8 +222,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
        if ((nrhs != 4))\r
                mexErrMsgTxt(\r
                                "erwarte (coordinates(Nx3), elements(Mx4), zeta(double), type(int))");\r
-       if (nlhs > 1)\r
-               mexErrMsgTxt("Funktion hat nur einen Rueckgabewert");\r
+       if (nlhs > 2 || nlhs == 0)\r
+               mexErrMsgTxt("Funktion hat maximal 2 Rueckgabewerte [A counts]");\r
 \r
        int cm = mxGetM(prhs[0]);\r
        int cn = mxGetN(prhs[0]);\r
@@ -264,9 +265,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
         * Variablen initialisieren\r
         */\r
        plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);\r
+\r
        double * A = mxGetPr(plhs[0]);\r
        double * C = mxGetPr(prhs[0]);\r
        double * E = mxGetPr(prhs[1]);\r
+       double * COUNT = NULL;\r
+\r
+       if(nlhs==2){\r
+               count = true;\r
+               plhs[1] = mxCreateDoubleMatrix(2,2,mxREAL);\r
+               COUNT = mxGetPr(plhs[1]);\r
+       }\r
 \r
        int type = (int) *(mxGetPr(prhs[3]));\r
        double * zeta = { 0 };\r
@@ -324,7 +333,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
         */\r
 #ifdef PARALLEL\r
        {\r
-#pragma omp parallel for schedule(static,1) private(i,j,k,tmp,x,y,xa,xb,ya,yb,d,dt,rx,rxa,rxb,ry,rya,ryb,lastRow,firstRow) shared(C,E,ctypeO,ctypeP,em,targetSize,actualNumberOfThreads)\r
+#pragma omp parallel for schedule(static,1) private(i,j,k,tmp,x,y,xa,xb,ya,yb,d,rx,rxa,rxb,ry,rya,ryb,lastRow,firstRow) shared(C,E,ctypeO,ctypeP,em,targetSize,actualNumberOfThreads)\r
 \r
                for (i = 0; i < actualNumberOfThreads; ++i) {\r
                        firstRow = (int) sqrt(2 * i * targetSize);\r
index 301260d67eb0d4848c141950a19fa7a138e12e54..8ee43f21cf06ccda449f8d51255a8d2d28d703b6 100644 (file)
@@ -367,21 +367,3 @@ function updateF2S(ele)
         G_ref_t(G_ref_f2sT) = 0;
     end
 end
-
-
-%% Berechnet die Flaecheninhalte der Elemente
-% function updateS(ele)
-% %Globale Variabelen aufbauen
-% global G_ref_E;
-% global G_ref_C;
-% global G_ref_s;
-%      for i = ele
-%     % normalized Vector on every triangle
-%     tri = G_ref_E(i,:);
-%     a = (G_ref_C(tri(2),:)-G_ref_C(tri(1),:));
-%     b = (G_ref_C(tri(4),:)-G_ref_C(tri(1),:));
-%     N = cross(a',b');
-%     N = N/norm(N);
-%     G_ref_s(i) = sqrt(sum(N.*N));
-%      end
-% end
\ No newline at end of file
index aa35c19429792c7a0a2870ba303c9da3c8f61ffb..c7e9b7fc847c2e0058671332009e3cb137d2d6e3 100644 (file)
 int GAUSS_size = GAUSS_SIZE[2];\r
 gauss * GAUSS_nodes = GAUSS_NODES[2];\r
 \r
+\r
+double count[][2] = {{0,0},{0,0}};\r
+\r
+\r
 /*\r
  * setzt den Quadraturgrad, also 2^q Auswertungsstellen\r
  */\r
@@ -29,6 +33,13 @@ int setQuad(int q) {
        return GAUSS_size;\r
 }\r
 \r
+/*\r
+ * gibt die count Matrix zurueck\r
+ */\r
+double ** getCount(){\r
+       return (double **) count;\r
+}\r
+\r
 /*\r
  * gibt das signum von x zurueck\r
  */\r
@@ -502,7 +513,8 @@ double cOrtO1(double b, double d, double t, double v, double d1, double d2,
 \r
 /*\r
  * fuer parallele Elemente\r
- * voll Analytisch\r
+ * volle Quadratur fuer zeta_Q zulaessige Elemente\r
+ * analytisch sonst\r
  */\r
 double cParO2(double b, double d, double t, double v, double d1, double d2,\r
                double d3, double* zeta) {\r
@@ -519,7 +531,7 @@ double cParO2(double b, double d, double t, double v, double d1, double d2,
 \r
 /*\r
  * fuer orthogonale Elemente\r
- * volle Quadratur fuer zeta_1 zulaessige Elemente\r
+ * volle Quadratur fuer zeta_Q zulaessige Elemente\r
  * analytisch sonst\r
  */\r
 double cOrtO2(double b, double d, double t, double v, double d1, double d2,\r
@@ -537,8 +549,8 @@ double cOrtO2(double b, double d, double t, double v, double d1, double d2,
 \r
 /*\r
  * fuer parallele Elemente\r
- * Quadratur ueber ein Element fuer zeta_2 zulaessige Elemente\r
- * volle Quadratur fuer zeta_1 zulaessige Elemente\r
+ * Quadratur ueber ein Element fuer zeta_E zulaessige Elemente\r
+ * volle Quadratur fuer zeta_Q zulaessige Elemente\r
  * analytisch sonst\r
  */\r
 double cParO3(double b, double d, double t, double v, double d1, double d2,\r
@@ -558,8 +570,8 @@ double cParO3(double b, double d, double t, double v, double d1, double d2,
 \r
 /*\r
  * fuer orthogonale Elemente\r
- * Quadratur ueber ein Element fuer zeta_2 zulaessige Elemente\r
- * volle Quadratur fuer zeta_1 zulaessige Elemente\r
+ * Quadratur ueber ein Element fuer zeta_E zulaessige Elemente\r
+ * volle Quadratur fuer zeta_Q zulaessige Elemente\r
  * analytisch sonst\r
  */\r
 double cOrtO3(double b, double d, double t, double v, double d1, double d2,\r
@@ -580,7 +592,7 @@ double cOrtO3(double b, double d, double t, double v, double d1, double d2,
 \r
 /*\r
  * fuer parallele Elemente\r
- * Quadratur ueber ein Element fuer zeta_2 zulaessige Elemente\r
+ * Quadratur ueber ein Element fuer zeta_E zulaessige Elemente\r
  * analytisch sonst\r
  */\r
 double cParO4(double b, double d, double t, double v, double d1, double d2,\r
@@ -589,14 +601,14 @@ double cParO4(double b, double d, double t, double v, double d1, double d2,
        //kurze Seite nach vorn\r
        switch_site_par(b, d, t, v, d1, d2);\r
 \r
-       if (zeta[0] * zeta[1] * (b * b + d * d) < dist2_par(b, d, t, v, d1, d2, d3))\r
+       if (zeta[1] * zeta[1] * (b * b + d * d) < dist2_par(b, d, t, v, d1, d2, d3))\r
                return intQ2A2<FY1Y2_par>(b, d, t, v, d1, d2, d3);\r
 \r
        return intA4<F_par>(b, d, t, v, d1, d2, d3);\r
 }\r
 /*\r
  * fuer orthogonale Elemente\r
- * Quadratur ueber ein Element fuer zeta_2 zulaessige Elemente\r
+ * Quadratur ueber ein Element fuer zeta_E zulaessige Elemente\r
  * analytisch sonst\r
  */\r
 double cOrtO4(double b, double d, double t, double v, double d1, double d2,\r
@@ -606,7 +618,7 @@ double cOrtO4(double b, double d, double t, double v, double d1, double d2,
        //kurze Seite nach vorn\r
        switch_site_ort(b, d, t, v, d1, d2, d3);\r
 \r
-       if (zeta[0] * zeta[1] * (b * b + d * d) < dist2_ort(b, d, t, v, d1, d2, d3))\r
+       if (zeta[1] * zeta[1] * (b * b + d * d) < dist2_ort(b, d, t, v, d1, d2, d3))\r
                return intQ2A2<FY2Y3_ort>(b, d, t, v, d1, d2, d3);\r
 \r
        return intA4<F_ort>(b, d, t, v, d1, d2, d3);\r
index 2dbf1e4ac9f3e4773ccdc292d6bbaec2b80a87b0..5a8516a57b87e792b1b7925ee12ba2a7b7b2a583 100644 (file)
@@ -20,22 +20,22 @@ double cParO1(double, double, double, double, double, double, double, double*);
 double cOrtO1(double, double, double, double, double, double, double, double*);
 
 /*
- * volle Quadratur fuer zeta_1 zulaessige Elemente
+ * volle Quadratur fuer zeta_Q zulaessige Elemente
  * analytisch sonst
  */
 double cParO2(double, double, double, double, double, double, double, double*);
 double cOrtO2(double, double, double, double, double, double, double, double*);
 
 /*
- * Quadratur ueber ein Element fuer zeta_2 zulaessige Elemente
- * volle Quadratur fuer zeta_1 zulaessige Elemente
+ * Quadratur ueber ein Element fuer zeta_E zulaessige Elemente
+ * volle Quadratur fuer zeta_Q zulaessige Elemente
  * analytisch sonst
  */
 double cParO3(double, double, double, double, double, double, double, double*);
 double cOrtO3(double, double, double, double, double, double, double, double*);
 
 /*
- * Quadratur ueber ein Element fuer zeta_1 zulaessige Elemente
+ * Quadratur ueber ein Element fuer zeta_E zulaessige Elemente
  * analytisch sonst
  */
 double cParO4(double, double, double, double, double, double, double, double*);
@@ -46,4 +46,9 @@ double cOrtO4(double, double, double, double, double, double, double, double*);
  */
 int setQuad(int);
 
+/*
+ * gibt die count Matrix zurueck
+ */
+double ** getCount();
+
 #endif
index 86193619d7ca41a2d64f1a38573a4f9ee4a18265..f867581553c16942d54b04415f0cb4438d12c585 100644 (file)
@@ -2,17 +2,18 @@
 format longG
 
 % Matrix MEX Funktion neu Compilieren
-mex -o 'mex_build_V.cpp' 'slpRectangle.cpp' CFLAGS="\$CFLAGS  -fopenmp" CXXFLAGS="\$CXXFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp"
+delete mex_build_V.mex*
+mex mex_build_V.cpp slpRectangle.cpp CFLAGS="\$CFLAGS -O2  -fopenmp" CXXFLAGS="\$CXXFLAGS -O2 -fopenmp" LDFLAGS="\$LDFLAGS -O2 -fopenmp"
 
 % mex mex_build_V.cpp slpRectangle.cpp
 
 % Test ausführen
 
 %Anzahl der Schritte oder wenn groeßer als 40 der Elemente
-steps = 11^4;
+steps = 4;
 
 %Art der Berechnungen
-type = [1 3 2];
+type = [1];
 zeta = { [2 2 2] [2 2 2] [2 2 2]};