From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Tue, 8 Nov 2011 10:27:10 +0000 (+0000)
Subject: [src] neue Tests
X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=e1770f91b544a864a339cd3c594040bd7c2f4d91;p=bacc.git

[src] neue Tests

git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@64 26120e32-c555-405d-b3e1-1f783fb42516
---

diff --git a/doc/doc1.pdf b/doc/doc1.pdf
index b0012de..3aa4b96 100644
Binary files a/doc/doc1.pdf and b/doc/doc1.pdf differ
diff --git a/doc/doc1.tex b/doc/doc1.tex
index c21c7dd..7b9937b 100644
--- a/doc/doc1.tex
+++ b/doc/doc1.tex
@@ -18,11 +18,15 @@
 
 \def\todo#1{\textcolor{red}{#1}}
 \def\Matlab{{\sc Matlab}}
-\def\q{\mathcal{Q}}
+\def\q{\Q}
 
 \def\Ea{$E$}
 \def\Eb{$\tilde E$}
 \def\v{\boldsymbol{v}}
+\def\D{\mathcal{D}}
+\def\Q{\mathcal{Q}}
+\def\G{\mathcal{G}}
+\def\L{\mathcal{L}}
 
 \newcommand{\dif}[2]{\frac{\partial#1}{\partial#2}}
 \newcommand{\enorm}[1]{\left| \! \left| \! \left|#1\right| \! \right| \! \right|}
@@ -94,9 +98,10 @@ Für $p \in \{1/2 , 0 , -1/2 , -1, -3/2\}$ gilt explizit:
 \begin{eqnarray*}
  G(p;y_1;y_2;x_1;x_2;\lambda) &:=& \int \int \{(x_1-y_1)^2 + (x_2-y_2)^2 + \lambda^2\}^p dy_2 dy_1
 \end{eqnarray*}
+Im Zuge der Berechnungen werden auch hier nur die Funktionen für bestimmte Parameter benötigt, die Untersuchung der restlichen werden wir hier nicht durchführen.
 
 \section{Integral über zwei Elemente}
-Bei der Integration über zwei Seitenelemente \Ea, \Eb~ muss man geometrisch zischen zwei Fällen unterscheiden. Entweder die beiden Elemente liegen in parallelen Ebenen oder in orthogonalen Ebenen.
+Bei der Integration über zwei Seitenelemente \Ea, \Eb~ haben wir geometrisch zischen zwei Fällen unterschieden. Entweder die beiden Elemente liegen in parallelen Ebenen oder in orthogonalen Ebenen.
 
 \subsection{Parallele Elemente}
 Liegen die beiden Elemente parallel zueinander lassen sie sich Folgendermaßen darstellen:
@@ -235,7 +240,7 @@ dy_2 dy_1 dx_2 dx_1\\
 \chapter{Semi-analytische Berechnung bestimmter Integrale}
 \section{Gauss-Quadratur}
 \begin{displaymath}
- \int_0^l f(x) dx \approx \sum_{k=0}^n w_kf(x_k) =: \mathcal{Q}(f)
+ \int_0^l f(x) dx \approx \sum_{k=0}^n w_kf(x_k) =: \Q(f)
 \end{displaymath}
 
 \section{Quadratur über ein Element}
@@ -245,15 +250,21 @@ dy_2 dy_1 dx_2 dx_1\\
 Das Integral über \Eb~ kann durch die analytische Doppelintegral Funktion ersetzt werden. Die Quadratur über \Ea~ wird nun mittels einer doppelten Gauss-Quadratur berechnet.
 \subsection{Zulässigkeitsbedingung}
 \begin{displaymath}
- \text{dist}(E,\tilde E) \req \mu min\{ \text{dist}
+  \D (E, \tilde E) \geq \mu \min\{ \G(E) , \G(\tilde E)\}
 \end{displaymath}
-
+Nur wenn die Zulässigkeitsbedingung gültig ist wird die Quadratur verwendet. Sonst wird Analytisch gerechnet.
 
 
 \section{Quadratur über eine Achse}
 \begin{displaymath}
  \int_{E} \int_{\tilde E} \approx \q_{I_1} \q_{J_1} \int_{I_2}\int_{J_2}
 \end{displaymath}
+\subsection{Zulässigkeitsbedingung}
+\begin{displaymath}
+  \D (I_1, J_1) \geq \mu \max\{ \L(I_2) , \L(J_2)\} 
+\end{displaymath}
+
+
 \section{Quadratur über eine Seite}
 \begin{displaymath}
  \int_{E} \int_{\tilde E} \approx \q_{I_1} \int_{I_2} \int_{J_1} \int_{J_2}
diff --git a/src/A_AnIso.m b/src/A_AnIso.m
index 56483be..d60bbcc 100644
--- a/src/A_AnIso.m
+++ b/src/A_AnIso.m
@@ -1,6 +1,6 @@
 function A_AnIso(file,times,mu,type,eta,eps)
-%A_ANISO Summary of this function goes here
-%   Detailed explanation goes here
+% A_AnIso(file,times,mu,type,eta,eps)
+
 
 A_loadMeshF(file)
 
diff --git a/src/A_Iso.m b/src/A_Iso.m
index e021601..bf51d6c 100644
--- a/src/A_Iso.m
+++ b/src/A_Iso.m
@@ -1,6 +1,6 @@
 function A_Iso(file,times,mu,type)
-%A_ANISO Summary of this function goes here
-%   DeTiiled explanation goes here
+%A_Iso(file,times,mu,type)
+
 
 A_loadMeshF(file)
 
diff --git a/src/A_stepAnIso.m b/src/A_stepAnIso.m
index beb3bd1..4b3fa1f 100644
--- a/src/A_stepAnIso.m
+++ b/src/A_stepAnIso.m
@@ -1,6 +1,5 @@
 function [dataAniso time er] = A_stepAnIso(mu,type,eta,eps)
-% dataAniso = A_stepAnIso(mu,type,eta,eps)
-%
+% [dataAniso time er] = A_stepAnIso(mu,type,eta,eps)
 
 global G_E;
 global G_C;
diff --git a/src/A_stepIso.m b/src/A_stepIso.m
index 0a7a22e..ad9bbe3 100644
--- a/src/A_stepIso.m
+++ b/src/A_stepIso.m
@@ -1,4 +1,5 @@
 function [dataIso timer] = A_stepIso(mu,type)
+% [dataIso timer] = A_stepIso(mu,type)
 
 global G_E;
 global G_C;
diff --git a/src/gauss.hpp b/src/gauss.hpp
new file mode 100644
index 0000000..17f9c1f
--- /dev/null
+++ b/src/gauss.hpp
@@ -0,0 +1,83 @@
+#ifndef GAUSS_NODES
+#define GAUSS_NODES
+
+#include <string.h>
+
+class gauss{
+public:
+	double * nod;
+	double * wei;
+};
+
+gauss* Quad = new gauss [6];
+
+void initeQuad(){
+
+double gn1[] = {0.2113248654051871,
+		0.7886751345948129
+		};
+
+double gw1[] = {0.4999999999999999,
+		0.4999999999999999
+		};
+
+
+
+Quad[0].nod = new double[1];
+Quad[0].nod[0] = 0.5;
+Quad[0].wei = new double[1];
+Quad[0].wei[0] = 1;
+
+Quad[1].nod = new double[2];
+memcpy(Quad[1].nod,gn1,sizeof(double)*2);
+Quad[1].wei = new double[2];
+memcpy(Quad[1].wei,gw1,sizeof(double)*2);
+
+Quad[2].nod = new double[4];
+memcpy(Quad[2].nod,gn4,sizeof(double)*4);
+Quad[2].wei = new double[4];
+memcpy(Quad[2].wei,gw4,sizeof(double)*4);
+
+}
+/*
+class Gauss {
+public:
+	Gauss(int n){
+		N = n;
+		node = new double[n];
+		weight = new double[n];
+
+
+	}
+	~Gauss(){
+		delete[] node;
+		delete[] weight;
+	}
+
+	int getN() const{
+		return N;
+	}
+
+	double getw(int i) const{
+		if(i<N)
+			return weight[i];
+		else
+			retrun 0;
+	}
+	double getn(int i) const {
+		if(i<N)
+			return node[i];
+		else
+			retrun 0;
+	}
+
+private:
+	int N;
+	double * node;
+	double * weight;
+};
+
+*/
+
+
+#endif
diff --git a/src/gauss.m b/src/gauss.m
index 795c59f..48a1099 100644
--- a/src/gauss.m
+++ b/src/gauss.m
@@ -30,7 +30,7 @@ beta = (1:n-1)./sqrt((2*(1:n-1)).^2-1);
 A = diag(beta,-1)+diag(beta,1);
 [eigenvector,nodes] = eig(A);
 [nodes,idx] = sort(diag(nodes));
-weights = 2*eigenvector(1,idx).^2; 
+weights = 2*eigenvector(1,idx).^2;
 
 if nargin >= 3
     a = varargin{1};
diff --git a/src/mex_build_AU.cpp b/src/mex_build_AU.cpp
index 3dc3aae..737851c 100644
--- a/src/mex_build_AU.cpp
+++ b/src/mex_build_AU.cpp
@@ -11,6 +11,8 @@
 #include "mex.h"
 #include <stdlib.h>
 
+//#include "gauss.hpp"
+
 #define M_EPS 10e-8
 
 //#include "tbb/parallel_for.h"
@@ -41,6 +43,9 @@ inline int dimOfThird(int a, int b) {
 
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
+//	initeQuad();
+//	cout << Quad[1].nod[0];
+
 	int i, j, k; //Schleifenindizes
 	double tmp; //Zwischenspeicherung der Einzelnen Werte
 	int count;
diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
index d984cf8..3dbe4c5 100644
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -6,6 +6,7 @@
 //#include <mex.h>
 
 #include "slpRectangle.hpp"
+//#include "gauss.hpp"
 
 using namespace std;
 
@@ -32,6 +33,10 @@ double inline max(double x, double y) {
 	return x < y ? y : x;
 }
 
+double inline min(double x, double y) {
+	return x > y ? y : x;
+}
+
 double g_QY(double p, double y, double x, double l) {
 
 	double sol = 0;
@@ -456,22 +461,6 @@ double calcParIntO3(double b, double d, double t, double v, double d1,
 double calcOrtIntO3(double b, double d, double t, double v, double d1,
 		double d2, double d3, double eta) {
 
-	//Elmente anordnen
-	if (d > t) {
-		double tmp = 0;
-
-		tmp = b;
-		b = v;
-		v = tmp;
-		tmp = d;
-		d = t;
-		t = tmp;
-
-		d1 = -d1;
-		d2 = -d2;
-		d3 = -d3;
-	}
-
 	if (max(b, t) * eta > d1) {
 		return 0;
 	} else {
diff --git a/src/test_calcInt1.m b/src/test_calcInt1.m
index 5a9e679..89515a5 100644
--- a/src/test_calcInt1.m
+++ b/src/test_calcInt1.m
@@ -1 +1 @@
-

elements=[1 2 3 4;5 6 7 8];
neigh = zeros(2,8);

dat = [];

h = 1;
coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 2 0 0; 2+h 0 0; 2+h h 0 ; 2 h 0];
figure(1)
current = coordinates(elements(1,[1:4,1])',:);
fill3(current(:,1),current(:,2),current(:,3),'g');
hold on
current = coordinates(elements(2,[1:4,1])',:);
fill3(current(:,1),current(:,2),current(:,3),'b');

h = h/2;
coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 2 0 0; 2+h 0 0; 2+h h 0 ; 2 h 0];
current = coordinates(elements(2,[1:4,1])',:);
fill3(current(:,1),current(:,2),current(:,3),'y');
hold off
legend('Element1(h/2)','Element2(h/2)','Element2(h/4)')

h = 2;
for I = 1:50

  h = h/2;
coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 2 0 0; 2+h 0 0; 2+h h 0 ; 2 h 0];
  
  A0 = mex_build_AU(coordinates,elements,0,0);
  A2 = mex_build_AU(coordinates,elements,1,2);
  A1 = mex_build_AU(coordinates,elements,1,1);
  I
  dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)];
    dat(dat<0)=0;
end

figure(2)
loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4)));

legend('Analytisch','Quad Element','Element vertauschen');

xlabel('Elementgroesse');
ylabel('Integral');
\ No newline at end of file
+

coo = @(h,diff)[0 0 0;1 0 0; 1 1 0;0 1  0; 0 0 0 ; 0+h 0 0; 0+h h 0 ; 0 h 0]+[zeros(4,3);repmat(diff,4,1)];
elements=[1 2 3 4;5 6 7 8];
neigh = zeros(2,8);
dat = [];

diff = [2 0 0];

%% Laage Übersicht
figure(1)
h = 1;
coordinates=coo(h,diff)
current = coordinates(elements(1,[1:4,1])',:);
fill3(current(:,1),current(:,2),current(:,3),'g');
hold on
current = coordinates(elements(2,[1:4,1])',:);
fill3(current(:,1),current(:,2),current(:,3),'b');

h = h/2;
coordinates=coo(h,diff);
current = coordinates(elements(2,[1:4,1])',:);
fill3(current(:,1),current(:,2),current(:,3),'y');
hold off
legend('Element1(h)','Element2(h)','Element2(h/2)');
title('Laage der Elemente');

%% Integrale bei kleiner werdenden Elementen

h = 2;
for I = 1:50

  h = h/2;
  coordinates=coo(h,diff);
  
  A0 = mex_build_AU(coordinates,elements,0,0);
  A2 = mex_build_AU(coordinates,elements,1,2);
  A1 = mex_build_AU(coordinates,elements,1,1);
  I
  dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)];
    dat(dat<0)=0;
end

figure(2)
loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4)));

legend('Analytisch','Quad Element','Element vertauschen','location','southeast');
xlabel('Elementgroesse (kürzeste Seite)');
ylabel('Integral');
title('Integral bei kleiner werdenden Element');

%% Netzverfeinerung ab einem Wert

h =  (1/2)^25;
coordinates=coo(h,diff);

A_loadMesh(coordinates,elements,neigh);
datA=[];
for I = 1:5
  datA(I,:) = A_stepIso(1,[0 1 3 2],1,0);
  I
end
figure(3)
loglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,6),datA(2:end,1),datA(2:end,12),datA(2:end,1),datA(2:end,9))

legend('Analytisch','Element vertauschen','Quad Element','Quad Achse','location','northeast');
xlabel('Anzahl der Elemente');
ylabel('mu Schaetzer');
title('mu Schaetzer mit "schlechtem" Startnetz');

figure(4)
loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,7),datA(:,1),datA(:,13),datA(:,1),datA(:,10))

legend('Analytisch','Element vertauschen','Quad Element','Quad Achse','location','southeast');
xlabel('Anzahl der Elemente');
ylabel('EnergieNorm ^2 ');
title('EnergieNorm ^2 mit "schlechtem" Startnetz');

datA
\ No newline at end of file
diff --git a/src/test_calcInt2.m b/src/test_calcInt2.m
index 9320aec..c86ccfc 100644
--- a/src/test_calcInt2.m
+++ b/src/test_calcInt2.m
@@ -5,8 +5,10 @@ neigh = zeros(2,8);
 
 dat = [];
 
+diff = 10;
+
 h = 1/2;
-coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 4 0 0; 4+h 0 0; 4+h 1 0 ; 4 1 0];
+coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; diff 0 0; diff+h 0 0; diff+h 1 0 ; diff 1 0];
 figure(1)
 current = coordinates(elements(1,[1:4,1])',:);
 fill3(current(:,1),current(:,2),current(:,3),'g');
@@ -15,7 +17,7 @@ current = coordinates(elements(2,[1:4,1])',:);
 fill3(current(:,1),current(:,2),current(:,3),'b');
 
 h = h/2;
-coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 4 0 0; 4+h 0 0; 4+h 1 0 ; 4 1 0];
+coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; diff 0 0; diff+h 0 0; diff+h 1 0 ; diff 1 0];
 current = coordinates(elements(2,[1:4,1])',:);
 fill3(current(:,1),current(:,2),current(:,3),'y');
 hold off
@@ -25,14 +27,14 @@ h = 2;
 for I = 1:50
 
   h = h/2;
-coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 4 0 0; 4+h 0 0; 4+h 1 0 ; 4 1 0];
+coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; diff 0 0; diff+h 0 0; diff+h 1 0 ; diff 1 0];
   
   A0 = mex_build_AU(coordinates,elements,0,0);
   A2 = mex_build_AU(coordinates,elements,1,2);
   A3 = mex_build_AU(coordinates,elements,1,3);
   I
   dat(I,1:4) = [h A0(1,2) A2(1,2) A3(1,2)];
-    dat(dat<0)=0;
+%    dat(dat<0)=0;
 end
 
 figure(2)
diff --git a/src/test_calcInt4.m b/src/test_calcInt4.m
index c2b60d7..7de2a56 100644
--- a/src/test_calcInt4.m
+++ b/src/test_calcInt4.m
@@ -33,7 +33,7 @@ coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 2 0 0; 2 0 h; 2 h h ; 2 h 0];
   I
   dat(I,1:4) = [h A0(1,2) A1(1,2) A2(1,2)];
   
-  dat(dat<0)=0;
+ % dat(dat<0)=0;
   
 end
 
diff --git a/src/vec2str.m b/src/vec2str.m
new file mode 100644
index 0000000..0f93a37
--- /dev/null
+++ b/src/vec2str.m
@@ -0,0 +1,14 @@
+function vec2str(vec)
+
+fprintf('{');
+
+for i = vec(1:end-1)
+  
+  fprintf('%6.16f,\n',i);
+  
+
+end
+
+fprintf('%6.16f\n',vec(end));
+
+fprintf('}\n');
\ No newline at end of file