\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|}
\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:
\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}
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}
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)
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)
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;
function [dataIso timer] = A_stepIso(mu,type)
+% [dataIso timer] = A_stepIso(mu,type)
global G_E;
global G_C;
--- /dev/null
+#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
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};
#include "mex.h"\r
#include <stdlib.h>\r
\r
+//#include "gauss.hpp"\r
+\r
#define M_EPS 10e-8\r
\r
//#include "tbb/parallel_for.h"\r
\r
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
\r
+// initeQuad();\r
+// cout << Quad[1].nod[0];\r
+\r
int i, j, k; //Schleifenindizes\r
double tmp; //Zwischenspeicherung der Einzelnen Werte\r
int count;\r
//#include <mex.h>\r
\r
#include "slpRectangle.hpp"\r
+//#include "gauss.hpp"\r
\r
using namespace std;\r
\r
return x < y ? y : x;\r
}\r
\r
+double inline min(double x, double y) {\r
+ return x > y ? y : x;\r
+}\r
+\r
double g_QY(double p, double y, double x, double l) {\r
\r
double sol = 0;\r
double calcOrtIntO3(double b, double d, double t, double v, double d1,\r
double d2, double d3, double eta) {\r
\r
- //Elmente anordnen\r
- if (d > t) {\r
- double tmp = 0;\r
-\r
- tmp = b;\r
- b = v;\r
- v = tmp;\r
- tmp = d;\r
- d = t;\r
- t = tmp;\r
-\r
- d1 = -d1;\r
- d2 = -d2;\r
- d3 = -d3;\r
- }\r
-\r
if (max(b, t) * eta > d1) {\r
return 0;\r
} else {\r
-\r\relements=[1 2 3 4;5 6 7 8];\rneigh = zeros(2,8);\r\rdat = [];\r\rh = 1;\rcoordinates=[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];\rfigure(1)\rcurrent = coordinates(elements(1,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'g');\rhold on\rcurrent = coordinates(elements(2,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'b');\r\rh = h/2;\rcoordinates=[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];\rcurrent = coordinates(elements(2,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'y');\rhold off\rlegend('Element1(h/2)','Element2(h/2)','Element2(h/4)')\r\rh = 2;\rfor I = 1:50\r\r h = h/2;\rcoordinates=[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];\r \r A0 = mex_build_AU(coordinates,elements,0,0);\r A2 = mex_build_AU(coordinates,elements,1,2);\r A1 = mex_build_AU(coordinates,elements,1,1);\r I\r dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)];\r dat(dat<0)=0;\rend\r\rfigure(2)\rloglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4)));\r\rlegend('Analytisch','Quad Element','Element vertauschen');\r\rxlabel('Elementgroesse');\rylabel('Integral');
\ No newline at end of file
+\r\rcoo = @(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)];\relements=[1 2 3 4;5 6 7 8];\rneigh = zeros(2,8);\rdat = [];\r\rdiff = [2 0 0];\r\r%% Laage Übersicht\rfigure(1)\rh = 1;\rcoordinates=coo(h,diff)\rcurrent = coordinates(elements(1,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'g');\rhold on\rcurrent = coordinates(elements(2,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'b');\r\rh = h/2;\rcoordinates=coo(h,diff);\rcurrent = coordinates(elements(2,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'y');\rhold off\rlegend('Element1(h)','Element2(h)','Element2(h/2)');\rtitle('Laage der Elemente');\r\r%% Integrale bei kleiner werdenden Elementen\r\rh = 2;\rfor I = 1:50\r\r h = h/2;\r coordinates=coo(h,diff);\r \r A0 = mex_build_AU(coordinates,elements,0,0);\r A2 = mex_build_AU(coordinates,elements,1,2);\r A1 = mex_build_AU(coordinates,elements,1,1);\r I\r dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)];\r dat(dat<0)=0;\rend\r\rfigure(2)\rloglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4)));\r\rlegend('Analytisch','Quad Element','Element vertauschen','location','southeast');\rxlabel('Elementgroesse (kürzeste Seite)');\rylabel('Integral');\rtitle('Integral bei kleiner werdenden Element');\r\r%% Netzverfeinerung ab einem Wert\r\rh = (1/2)^25;\rcoordinates=coo(h,diff);\r\rA_loadMesh(coordinates,elements,neigh);\rdatA=[];\rfor I = 1:5\r datA(I,:) = A_stepIso(1,[0 1 3 2],1,0);\r I\rend\rfigure(3)\rloglog(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))\r\rlegend('Analytisch','Element vertauschen','Quad Element','Quad Achse','location','northeast');\rxlabel('Anzahl der Elemente');\rylabel('mu Schaetzer');\rtitle('mu Schaetzer mit "schlechtem" Startnetz');\r\rfigure(4)\rloglog(datA(:,1),datA(:,4),datA(:,1),datA(:,7),datA(:,1),datA(:,13),datA(:,1),datA(:,10))\r\rlegend('Analytisch','Element vertauschen','Quad Element','Quad Achse','location','southeast');\rxlabel('Anzahl der Elemente');\rylabel('EnergieNorm ^2 ');\rtitle('EnergieNorm ^2 mit "schlechtem" Startnetz');\r\rdatA\r
\ No newline at end of file
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');
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
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)
I
dat(I,1:4) = [h A0(1,2) A1(1,2) A2(1,2)];
- dat(dat<0)=0;
+ % dat(dat<0)=0;
end
--- /dev/null
+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
projects / bacc.git / commitdiff