</toolChain>
</folderInfo>
<sourceEntries>
- <entry excluding="src/bem3d" flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="sourcePath" name=""/>
+ <entry excluding="src/maik|src/bem3d" flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="sourcePath" name=""/>
</sourceEntries>
</configuration>
</storageModule>
</toolChain>
</folderInfo>
<sourceEntries>
- <entry excluding="src/bem3d" flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="sourcePath" name=""/>
+ <entry excluding="src/maik|src/bem3d" flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="sourcePath" name=""/>
</sourceEntries>
</configuration>
</storageModule>
--- /dev/null
+\documentclass[a4paper,10pt,fleqn]{report}
+
+\usepackage{fullpage} %Seiten etwas Größer
+\usepackage{amsmath} %Mathematische Symbole
+%\usepackage{moreverb}
+\usepackage{graphicx,psfrag,subfig} %Grafiken einbinden/Texte ersetzen/Bilder nebeneinander
+%\usepackage{ifthen}
+%\usepackage{showkeys}
+%\usepackage[final,numbered]{mcode}
+\usepackage{color} %Farben benutzen und Definieren
+\usepackage{hyperref} %Links im Inhaltsverzeichnis
+\hypersetup{linkbordercolor={1 1 1},citebordercolor={1 1 1},urlbordercolor={1 1 1}}
+
+
+\usepackage[ngerman]{babel} %Sprachpacket für Überschriften
+\usepackage[utf8]{inputenc} %Eingabekodierung
+\usepackage{fixltx2e} %Deutschsprach Bugs
+
+
+\def\todo#1{\textcolor{red}{#1}}
+\def\Matlab{{\sc Matlab}}
+
+
+\date{29.08.2011}
+\author{P. Schaefer}
+
+\begin{document}
+\tableofcontents
+\clearpage
+\chapter{Einleitung}
+\section{Problem}
+Zwei Elemente $E = v + [(0,l_1) \times (0,l_2) \times \{0\}]$ und $\hat E = \hat v + [(0,\hat l_1) \times (0,\hat l_2) \times \{0\}]$
+
+\noindent
+Berechnet werden soll:
+\begin{eqnarray*}
+ -\frac{1}{4\pi} \int_E \int_{\hat E} \frac{1}{|x-y|} ds_y ds_x
+\end{eqnarray*}
+
+\section{Gitter}
+Alle Gitter bestehen nur aus achsenorientierten Rechtecken.
+Entscheidend für die Tests sind Gitter in einer Ebene und Ganze Körper...
+
+\chapter{Analytisch}
+
+\section{einfach Inetgral}
+\begin{eqnarray*}
+ g(p;y;x;\lambda) &:=& \int \{(x-y)^2 + \lambda^2 \}^p dy
+\end{eqnarray*}
+
+
+\section{doppel Integral}
+\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*}
+
+
+
+\section{Parallele Elemente}
+Für parallele Elemente lässt sich zeigen:
+\begin{eqnarray*}
+F(x_1,x_2,y_1,y_2,\delta_1,\delta_2,\delta_3) &:=& \int \int \int \int \left\{(x_1-y_1-\delta_1)^2+(x_2-y_2-\delta_2)^2+\delta_3^2\right\}^{-1/2}dy_2 dy_1 dx_2 dx_1
+\end{eqnarray*}
+
+\section{Ortogonale Elemente}
+Für ortogonale Elemente lässt sich zeigen:
+\begin{eqnarray*}
+F(x_1,x_2,y_1,y_2,\delta_1,\delta_2,\delta_3) &:=& \int \int \int \int \left\{(x_1-\delta_1)^2+(x_2-y_2-\delta_2)^2+(-y_3-\delta_3)^2\right\}^{-1/2}dy_3 dy_2 dx_2 dx_1
+\end{eqnarray*}
+
+\chapter{Quadratur}
+\section{Was ist das}
+gauss-Quadratur
+\section{Umsetzung für BLA}
+
+
+
+\chapter{Schätzer und Kreterien}
+\section{mögliche Laagen der Elemente}
+so so oder so...
+\section{Fehlerschätzer}
+\section{möglichkeiten für Kreterien}
+
+
+\end{document}
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
- plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
-// plhs[0] = mxCreateDoubleMatrix(1, 2, mxREAL);
-
-
- //printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-
- //sol = G00(p,y1,y2,x1,x2,l);
*mxGetPr(plhs[0]) = F_ort(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
- *mxGetPr(plhs[1]) = FLO_perp(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
-// printf("%f",FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6])));
-
-// *(mxGetPr(plhs[0])+1) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
}
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
- plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
-// plhs[0] = mxCreateDoubleMatrix(1, 2, mxREAL);
-
-
- //printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-
- //sol = G00(p,y1,y2,x1,x2,l);
*mxGetPr(plhs[0]) = F_par(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
- *mxGetPr(plhs[1]) = F_parY1Y2_2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
-// printf("%f",FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6])));
-
-// *(mxGetPr(plhs[0])+1) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
}
+++ /dev/null
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.hpp"
-
-#define my_PI 3.141592653589793
-#define EPS 0.02
-
-using namespace std;
-
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
- //sicherheitsabfragen zu Datengroessen
- if (nrhs != 7)
- mexErrMsgTxt("expected (x1,x2,y1,y2,d1,d2,d3)");
- if (nlhs > 2)
- mexErrMsgTxt("has maximal two output arguments");
-
-
- plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
-// plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
-// plhs[0] = mxCreateDoubleMatrix(1, 2, mxREAL);
-
-
- //printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-
- //sol = G00(p,y1,y2,x1,x2,l);
-
- *mxGetPr(plhs[0]) = F_parQY2X2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-// *mxGetPr(plhs[1]) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
-// printf("%f",FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6])));
-
-// *(mxGetPr(plhs[0])+1) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
-
-}
-
-
-
+++ /dev/null
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.hpp"
-
-#define my_PI 3.141592653589793
-#define EPS 0.02
-
-using namespace std;
-
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
- //sicherheitsabfragen zu Datengroessen
- if (nrhs != 7)
- mexErrMsgTxt("expected (x1,x2,y1,y2,d1,d2,d3)");
- if (nlhs > 2)
- mexErrMsgTxt("has maximal two output arguments");
-
-
- plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
-// plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
-// plhs[0] = mxCreateDoubleMatrix(1, 2, mxREAL);
-
-
- //printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-
- //sol = G00(p,y1,y2,x1,x2,l);
-
- *mxGetPr(plhs[0]) = F_parY1Y2_2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-// *mxGetPr(plhs[1]) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
-// printf("%f",FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6])));
-
-// *(mxGetPr(plhs[0])+1) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
-
-}
-
-
-
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
- //printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-
- //sol = G00(p,y1,y2,x1,x2,l);
-
- *mxGetPr(plhs[0]) = G00(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
- *mxGetPr(plhs[1]) = doubleQuad(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
-
+ *mxGetPr(plhs[0]) = G_AY1Y2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
+ *mxGetPr(plhs[1]) = G_QY1Y2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
}
// y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
tmp\r
- = apply0Int4(F_par, fabs(xa[rxa]), fabs(xb[rxb]),\r
+ = calcParIntA( fabs(xa[rxa]), fabs(xb[rxb]),\r
fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],\r
d[rxb], d[rx]);\r
// printf("%d%d|%.2f\n",j,k,tmp);\r
// y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
tmp\r
- = apply0Int4(F_par, fabs(xa[rxa]), fabs(xb[rxb]),\r
+ = calcParIntA( fabs(xa[rxa]), fabs(xb[rxb]),\r
fabs(yb[rxa]), fabs(ya[rxb]), d[rxa],\r
d[rxb], d[rx]);\r
// printf("%d%d|%.2f\n",j,k,tmp);\r
\r
if (rxa == rya) {\r
tmp\r
- = apply0Int4(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),\r
+ = calcOrtIntA( fabs(xb[rxb]), fabs(xa[rxa]),\r
fabs(ya[rya]), fabs(yb[ryb]), d[rxb],\r
d[rxa], d[rx]);\r
} else if (rxa == ryb) {\r
tmp\r
- = apply0Int4(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),\r
+ = calcOrtIntA( fabs(xb[rxb]), fabs(xa[rxa]),\r
fabs(yb[ryb]), fabs(ya[rya]), d[rxb],\r
d[rxa], d[rx]);\r
} else if (rxb == rya) {\r
tmp\r
- = apply0Int4(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),\r
+ = calcOrtIntA( fabs(xa[rxa]), fabs(xb[rxb]),\r
fabs(ya[rya]), fabs(yb[ryb]), d[rxa],\r
d[rxb], d[rx]);\r
} else {\r
tmp\r
- = apply0Int4(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),\r
+ = calcOrtIntA( fabs(xa[rxa]), fabs(xb[rxb]),\r
fabs(yb[ryb]), fabs(ya[rya]), d[rxa],\r
d[rxb], d[rx]);\r
}\r
- A[(k * em) + j] = tmp;\r
\r
}\r
A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
\r
//LageInformationen\r
int rx, rxa, rxb, ry, rya, ryb;\r
- int rtxa, rtxb, rtx;\r
\r
//Ausrechnen\r
for (j = 0; j < em; ++j) {\r
// Lageeigenschaften des Flaechenstuecks\r
//if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
\r
- rtx = dimOfVec(xn);\r
+ rx = dimOfVec(xn);\r
\r
- rtxa = dimOfVec(xa);\r
+ rxa = dimOfVec(xa);\r
\r
- rtxb = dimOfVec(xb);\r
+ rxb = dimOfVec(xb);\r
\r
\r
//kleinste Ecke finden und fuer \delta verwenden\r
\r
- if (xa[rtxa] > 0) {\r
- if (xb[rtxb] > 0) {\r
+ if (xa[rxa] > 0) {\r
+ if (xb[rxb] > 0) {\r
for (i = 0; i < 3; ++i)\r
dt[i] = -x0[i];\r
} else {\r
dt[i] = -x3[i];\r
}\r
} else {\r
- if (xb[rtxb] > 0) {\r
+ if (xb[rxb] > 0) {\r
for (i = 0; i < 3; ++i)\r
dt[i] = -x1[i];\r
} else {\r
}\r
}\r
\r
- //for (i=0;i<3;++i)\r
- // dt[i] = 0;\r
-\r
-\r
- // printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);\r
+ xda = xa[rxa];\r
+ xdb = xb[rxb];\r
\r
for (k = 0; k < em; ++k) {\r
\r
- xda = xa[rtxa];\r
- xdb = xb[rtxb];\r
-\r
- rx = rtx;\r
- rxa = rtxa;\r
- rxb = rtxb;\r
\r
\r
y0[0] = C[(int) E[k] - 1];\r
yda = ya[rya];\r
ydb = yb[ryb];\r
\r
-\r
- if(xda*xda+xdb*xdb>yda*yda+ydb*ydb){\r
-\r
- /* printf("+(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
- printf("[%.2f,%.2f@%.2f,%.2f]",xda,xb[rxb],ya[rya],yb[ryb]);\r
- printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);\r
- printf("\n");*/\r
-\r
-\r
- dtmp = xda; xda = yda; yda = dtmp;\r
- dtmp = xdb; xdb = ydb; ydb = dtmp;\r
-\r
- itmp = rxa; rxa = rya; rya = itmp;\r
- itmp = rxb; rxb = ryb; ryb = itmp;\r
- itmp = rx; rx = ry; ry = itmp;\r
-\r
- for (i = 0; i<3;++i)\r
- d[i] = -d[i];\r
-\r
-\r
-\r
- /* printf("-(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
- printf("[%.2f,%.2f@%.2f,%.2f]",xda,xb[rxb],ya[rya],yb[ryb]);\r
- printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);\r
- printf("\n");*/\r
- }\r
-\r
-\r
//printf("(%d,%d)",rx,ry);\r
if (rx == ry) {\r
//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
// y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
\r
- // if(xda*xda+xdb*xdb<=yda*yda+ydb*ydb){\r
-\r
tmp\r
- = apply0Int4(F_par, fabs(xda), fabs(xdb),\r
+ = calcParIntO1( fabs(xda), fabs(xdb),\r
fabs(yda), fabs(ydb), d[rxa],\r
d[rxb], d[rx]);\r
\r
- // }else{\r
-\r
- // tmp\r
- // = apply0Int4(F_par, fabs(yda), fabs(ydb),\r
- // fabs(xda), fabs(xdb), -d[rxa],\r
- // -d[rxb], -d[rx]);\r
- // }\r
- // if(fabs(tmp-tmp2)>10e-16)\r
- // printf("wtf");\r
- // printf("%d%d|%.2f\n",j,k,tmp);\r
} else {\r
//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
// y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
tmp\r
- = apply0Int4(F_par, fabs(xda), fabs(xdb),\r
+ = calcParIntO1( fabs(xda), fabs(xdb),\r
fabs(ydb), fabs(yda), d[rxa],\r
d[rxb], d[rx]);\r
// printf("%d%d|%.2f\n",j,k,tmp);\r
\r
if (rxa == rya) {\r
tmp\r
- = apply0Int4(F_ort, fabs(xdb), fabs(xda),\r
+ = calcOrtIntO1( fabs(xdb), fabs(xda),\r
fabs(yda), fabs(ydb), d[rxb],\r
d[rxa], d[rx]);\r
} else if (rxa == ryb) {\r
tmp\r
- = apply0Int4(F_ort, fabs(xdb), fabs(xda),\r
+ = calcOrtIntO1( fabs(xdb), fabs(xda),\r
fabs(ydb), fabs(yda), d[rxb],\r
d[rxa], d[rx]);\r
} else if (rxb == rya) {\r
tmp\r
- = apply0Int4(F_ort, fabs(xda), fabs(xdb),\r
+ = calcOrtIntO1( fabs(xda), fabs(xdb),\r
fabs(yda), fabs(ydb), d[rxa],\r
d[rxb], d[rx]);\r
} else {\r
tmp\r
- = apply0Int4(F_ort, fabs(xda), fabs(xdb),\r
+ = calcOrtIntO1( fabs(xda), fabs(xdb),\r
fabs(ydb), fabs(yda), d[rxa],\r
d[rxb], d[rx]);\r
}\r
/*\r
* voll Analytisch\r
* Element mit groesserer Diagonale nach aussen gedreht\r
- *\r
+ * Quadratur ueber E_j bei dist(E_j,E_k) >= eta*dia(E_j) | dia(E_j) <= dia(E_k)\r
*\r
*/\r
\r
\r
#include "slpRectangle.hpp"\r
\r
-#define EPS 0.02\r
+//Quadraturverwendung!!!\r
+#define eta 0.8\r
\r
using namespace std;\r
\r
\r
//LageInformationen\r
int rx, rxa, rxb, ry, rya, ryb;\r
- int rtxa, rtxb, rtx;\r
\r
//Ausrechnen\r
for (j = 0; j < em; ++j) {\r
// Lageeigenschaften des Flaechenstuecks\r
//if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
\r
- rtx = dimOfVec(xn);\r
+ rx = dimOfVec(xn);\r
\r
- rtxa = dimOfVec(xa);\r
+ rxa = dimOfVec(xa);\r
\r
- rtxb = dimOfVec(xb);\r
+ rxb = dimOfVec(xb);\r
\r
\r
//kleinste Ecke finden und fuer \delta verwenden\r
\r
- if (xa[rtxa] > 0) {\r
- if (xb[rtxb] > 0) {\r
+ if (xa[rxa] > 0) {\r
+ if (xb[rxb] > 0) {\r
for (i = 0; i < 3; ++i)\r
dt[i] = -x0[i];\r
} else {\r
dt[i] = -x3[i];\r
}\r
} else {\r
- if (xb[rtxb] > 0) {\r
+ if (xb[rxb] > 0) {\r
for (i = 0; i < 3; ++i)\r
dt[i] = -x1[i];\r
} else {\r
}\r
}\r
\r
- //for (i=0;i<3;++i)\r
- // dt[i] = 0;\r
-\r
-\r
- // printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);\r
+ xda = xa[rxa];\r
+ xdb = xb[rxb];\r
\r
for (k = 0; k < em; ++k) {\r
\r
- xda = xa[rtxa];\r
- xdb = xb[rtxb];\r
-\r
- rx = rtx;\r
- rxa = rtxa;\r
- rxb = rtxb;\r
\r
\r
y0[0] = C[(int) E[k] - 1];\r
yda = ya[rya];\r
ydb = yb[ryb];\r
\r
- //groesseres Element nach aussen drehen\r
- if(xda*xda+xdb*xdb>yda*yda+ydb*ydb){\r
-\r
- dtmp = xda; xda = yda; yda = dtmp;\r
- dtmp = xdb; xdb = ydb; ydb = dtmp;\r
-\r
- itmp = rxa; rxa = rya; rya = itmp;\r
- itmp = rxb; rxb = ryb; ryb = itmp;\r
- itmp = rx; rx = ry; ry = itmp;\r
-\r
- for (i = 0; i<3;++i)\r
- d[i] = -d[i];\r
-\r
- }\r
-\r
- //if(xda*xda+xdb*xdb>=d[0]*d[0]+d[1]*d[1]+d[2]*d[2]){\r
- // printf(".");\r
- //}\r
-\r
//printf("(%d,%d)",rx,ry);\r
if (rx == ry) {\r
//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
// y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
\r
- // if(xda*xda+xdb*xdb<d[0]*d[0]+d[1]*d[1]+d[2]*d[2]){\r
-\r
tmp\r
- = apply0Int4(F_par, fabs(xda), fabs(xdb),\r
+ = calcParIntO2( fabs(xda), fabs(xdb),\r
fabs(yda), fabs(ydb), d[rxa],\r
- d[rxb], d[rx]);\r
+ d[rxb], d[rx],eta);\r
\r
- // }else{\r
-\r
- tmp2\r
- = apply0Int4(F_parY1Y2_2, fabs(xda), fabs(xdb),\r
- fabs(yda), fabs(ydb), d[rxa],\r
- d[rxb], d[rx]);\r
- // }\r
- if(fabs(tmp-tmp2)>10e-8)\r
- return;\r
- // printf("%d%d|%.2f\n",j,k,tmp);\r
} else {\r
//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
// y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
tmp\r
- = apply0Int4(F_par, fabs(xda), fabs(xdb),\r
+ = calcParIntO2( fabs(xda), fabs(xdb),\r
fabs(ydb), fabs(yda), d[rxa],\r
- d[rxb], d[rx]);\r
+ d[rxb], d[rx],eta);\r
// printf("%d%d|%.2f\n",j,k,tmp);\r
}\r
\r
\r
if (rxa == rya) {\r
tmp\r
- = apply0Int4(F_ort, fabs(xdb), fabs(xda),\r
+ = calcOrtIntO2( fabs(xdb), fabs(xda),\r
fabs(yda), fabs(ydb), d[rxb],\r
- d[rxa], d[rx]);\r
+ d[rxa], d[rx],eta);\r
} else if (rxa == ryb) {\r
tmp\r
- = apply0Int4(F_ort, fabs(xdb), fabs(xda),\r
+ = calcOrtIntO2( fabs(xdb), fabs(xda),\r
fabs(ydb), fabs(yda), d[rxb],\r
- d[rxa], d[rx]);\r
+ d[rxa], d[rx],eta);\r
} else if (rxb == rya) {\r
tmp\r
- = apply0Int4(F_ort, fabs(xda), fabs(xdb),\r
+ = calcOrtIntO2( fabs(xda), fabs(xdb),\r
fabs(yda), fabs(ydb), d[rxa],\r
- d[rxb], d[rx]);\r
+ d[rxb], d[rx],eta);\r
} else {\r
tmp\r
- = apply0Int4(F_ort, fabs(xda), fabs(xdb),\r
+ = calcOrtIntO2( fabs(xda), fabs(xdb),\r
fabs(ydb), fabs(yda), d[rxa],\r
- d[rxb], d[rx]);\r
+ d[rxb], d[rx],eta);\r
}\r
\r
}\r
--- /dev/null
+/*\r
+ * voll Analytisch\r
+ * Element mit groesserer Diagonale nach aussen gedreht\r
+ * Quadratur ueber E_j bei dist(E_j,E_k) >= eta*dia(E_j) | dia(E_j) <= dia(E_k)\r
+ *\r
+ */\r
+\r
+\r
+#include <iostream>\r
+#include <cmath>\r
+#include <cassert>\r
+#include "mex.h"\r
+#include <stdlib.h>\r
+\r
+//#include "tbb/parallel_for.h"\r
+\r
+\r
+#include "slpRectangle.hpp"\r
+\r
+//Quadraturverwendung!!!\r
+#define eta 0.8\r
+\r
+using namespace std;\r
+\r
+inline int dimOfVec(double * vec)\r
+{\r
+ if (vec[2] != 0)\r
+ return 2;\r
+ else if (vec[1] != 0)\r
+ return 1;\r
+ else\r
+ return 0;\r
+}\r
+\r
+\r
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
+\r
+ int i, j, k;\r
+\r
+ //sicherheitsabfragen zu Datengroessen\r
+ if (nrhs != 2)\r
+ mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx4))");\r
+ if (nlhs > 1)\r
+ mexErrMsgTxt("has only one output argument");\r
+\r
+ int cm = mxGetM(prhs[0]);\r
+ int cn = mxGetN(prhs[0]);\r
+\r
+ if (cn != 3)\r
+ mexErrMsgTxt("expected coordinates (Nx3)");\r
+ int em = mxGetM(prhs[1]);\r
+ int en = mxGetN(prhs[1]);\r
+ if (en != 4)\r
+ mexErrMsgTxt("expected elements (Mx4)");\r
+ //Vorbereitung der Daten\r
+\r
+ plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);\r
+ double * A = mxGetPr(plhs[0]);\r
+ double * C = mxGetPr(prhs[0]);\r
+ double * E = mxGetPr(prhs[1]);\r
+\r
+ double * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * xn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+\r
+ double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * yn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+\r
+ double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+\r
+ double tmp,tmp2;\r
+\r
+ int itmp;\r
+ double dtmp;\r
+ double * vtmp;\r
+\r
+ long double xda,xdb,yda,ydb;\r
+\r
+ //LageInformationen\r
+ int rx, rxa, rxb, ry, rya, ryb;\r
+\r
+ //Ausrechnen\r
+ for (j = 0; j < em; ++j) {\r
+ x0[0] = C[(int) E[j] - 1];\r
+ x0[1] = C[cm + (int) E[j] - 1];\r
+ x0[2] = C[2 * cm + (int) E[j] - 1];\r
+\r
+ x1[0] = C[(int) E[em + j] - 1];\r
+ x1[1] = C[cm + (int) E[em + j] - 1];\r
+ x1[2] = C[2 * cm + (int) E[em + j] - 1];\r
+\r
+ x2[0] = C[(int) E[2 * em + j] - 1];\r
+ x2[1] = C[cm + (int) E[2 * em + j] - 1];\r
+ x2[2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
+\r
+ x3[0] = C[(int) E[3 * em + j] - 1];\r
+ x3[1] = C[cm + (int) E[3 * em + j] - 1];\r
+ x3[2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
+\r
+ for (i = 0; i < 3; ++i)\r
+ xa[i] = x1[i] - x0[i];\r
+\r
+ for (i = 0; i < 3; ++i)\r
+ xb[i] = x3[i] - x0[i];\r
+\r
+ //Kreuzprodukt\r
+ xn[0] = (xa[1]) * (xb[2]) - (xa[2]) * (xb[1]);\r
+ xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);\r
+ xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);\r
+\r
+ // Lageeigenschaften des Flaechenstuecks\r
+ //if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
+\r
+ rx = dimOfVec(xn);\r
+\r
+ rxa = dimOfVec(xa);\r
+\r
+ rxb = dimOfVec(xb);\r
+\r
+\r
+ //kleinste Ecke finden und fuer \delta verwenden\r
+\r
+ if (xa[rxa] > 0) {\r
+ if (xb[rxb] > 0) {\r
+ for (i = 0; i < 3; ++i)\r
+ dt[i] = -x0[i];\r
+ } else {\r
+ for (i = 0; i < 3; ++i)\r
+ dt[i] = -x3[i];\r
+ }\r
+ } else {\r
+ if (xb[rxb] > 0) {\r
+ for (i = 0; i < 3; ++i)\r
+ dt[i] = -x1[i];\r
+ } else {\r
+ for (i = 0; i < 3; ++i)\r
+ dt[i] = -x2[i];\r
+ //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
+ }\r
+ }\r
+\r
+ xda = xa[rxa];\r
+ xdb = xb[rxb];\r
+\r
+ for (k = 0; k < em; ++k) {\r
+\r
+\r
+\r
+ y0[0] = C[(int) E[k] - 1];\r
+ y0[1] = C[cm + (int) E[k] - 1];\r
+ y0[2] = C[2 * cm + (int) E[k] - 1];\r
+\r
+ y1[0] = C[(int) E[em + k] - 1];\r
+ y1[1] = C[cm + (int) E[em + k] - 1];\r
+ y1[2] = C[2 * cm + (int) E[em + k] - 1];\r
+\r
+ y2[0] = C[(int) E[2 * em + k] - 1];\r
+ y2[1] = C[cm + (int) E[2 * em + k] - 1];\r
+ y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
+\r
+ y3[0] = C[(int) E[3 * em + k] - 1];\r
+ y3[1] = C[cm + (int) E[3 * em + k] - 1];\r
+ y3[2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
+\r
+ for (i = 0; i < 3; ++i)\r
+ ya[i] = y1[i] - y0[i];\r
+\r
+ for (i = 0; i < 3; ++i)\r
+ yb[i] = y3[i] - y0[i];\r
+\r
+ yn[0] = (ya[1]) * (yb[2]) - (ya[2]) * (yb[1]);\r
+ yn[1] = (ya[2]) * (yb[0]) - (ya[0]) * (yb[2]);\r
+ yn[2] = (ya[0]) * (yb[1]) - (ya[1]) * (yb[0]);\r
+\r
+ ry = dimOfVec(yn);\r
+\r
+ rya = dimOfVec(ya);\r
+\r
+ ryb = dimOfVec(yb);\r
+\r
+ //kleinste Ecke finden und fuer \delta verwenden\r
+\r
+ if (ya[rya] > 0) {\r
+ if (yb[ryb] > 0) {\r
+ for (i = 0; i < 3; ++i)\r
+ d[i] = dt[i] + y0[i];\r
+ } else {\r
+ for (i = 0; i < 3; ++i)\r
+ d[i] = dt[i] + y3[i];\r
+ }\r
+ } else {\r
+ if (yb[ryb] > 0) {\r
+ for (i = 0; i < 3; ++i)\r
+ d[i] = dt[i] + y1[i];\r
+ } else {\r
+ for (i = 0; i < 3; ++i)\r
+ d[i] = dt[i] + y2[i];\r
+ //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
+ }\r
+ }\r
+\r
+ yda = ya[rya];\r
+ ydb = yb[ryb];\r
+\r
+ //printf("(%d,%d)",rx,ry);\r
+ if (rx == ry) {\r
+ //printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
+ //printf("[%.2f,%.2f@%.2f,%.2f]",xda,xdb,yda,ydb);\r
+ //printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);\r
+ //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
+ if (rxa == rya) {\r
+ //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
+ // y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
+ //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
+\r
+ tmp\r
+ = calcParIntO3( fabs(xda), fabs(xdb),\r
+ fabs(yda), fabs(ydb), d[rxa],\r
+ d[rxb], d[rx],eta);\r
+\r
+ } else {\r
+ //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
+ // y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
+ //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
+ tmp\r
+ = calcParIntO3( fabs(xda), fabs(xdb),\r
+ fabs(ydb), fabs(yda), d[rxa],\r
+ d[rxb], d[rx],eta);\r
+ // printf("%d%d|%.2f\n",j,k,tmp);\r
+ }\r
+\r
+ } else {\r
+ //printf("(%d,%d@%d,%d)", rxa, rxb, rya, ryb);\r
+ //printf("(%d,%d)", rx, ry);\r
+ //printf("\n");\r
+\r
+ if (rxa == rya) {\r
+ tmp\r
+ = calcOrtIntO3( fabs(xdb), fabs(xda),\r
+ fabs(yda), fabs(ydb), d[rxb],\r
+ d[rxa], d[rx],eta);\r
+ } else if (rxa == ryb) {\r
+ tmp\r
+ = calcOrtIntO3( fabs(xdb), fabs(xda),\r
+ fabs(ydb), fabs(yda), d[rxb],\r
+ d[rxa], d[rx],eta);\r
+ } else if (rxb == rya) {\r
+ tmp\r
+ = calcOrtIntO3( fabs(xda), fabs(xdb),\r
+ fabs(yda), fabs(ydb), d[rxa],\r
+ d[rxb], d[rx],eta);\r
+ } else {\r
+ tmp\r
+ = calcOrtIntO3( fabs(xda), fabs(xdb),\r
+ fabs(ydb), fabs(yda), d[rxa],\r
+ d[rxb], d[rx],eta);\r
+ }\r
+\r
+ }\r
+ A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
+\r
+ }\r
+\r
+ }\r
+ //printf("\n");\r
+ //Rueckgabe (eventuell zurueck schreiben)\r
+ // mxFree(x0);\r
+ //mxFree(x1);\r
+ //mxFree(x3);\r
+ //mxFree(y0);\r
+ //mxFree(y1);\r
+ //mxFree(y3);\r
+ //mxFree(xn);\r
+ //mxFree(yn);\r
+ //mxFree(d);\r
+\r
+ return;\r
+}\r
--- /dev/null
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+#include "slpRectangle.hpp"
+
+#define my_PI 3.141592653589793
+#define EPS 0.02
+
+using namespace std;
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+ //sicherheitsabfragen zu Datengroessen
+ if (nrhs != 7)
+ mexErrMsgTxt("expected (x1,x2,y1,y2,d1,d2,d3)");
+ if (nlhs > 2)
+ mexErrMsgTxt("has maximal two output arguments");
+
+ plhs[0] = mxCreateDoubleMatrix(1, 3, mxREAL);
+
+ *mxGetPr(plhs[0]) = calcortIntA(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+// *(mxGetPr(plhs[0])+1) = calcParIntQX1X2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+// *(mxGetPr(plhs[0])+2) = calcParIntQY1X1(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+
+}
+
+
+
--- /dev/null
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+#include "slpRectangle.hpp"
+
+#define my_PI 3.141592653589793
+#define EPS 0.02
+
+using namespace std;
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+ //sicherheitsabfragen zu Datengroessen
+ if (nrhs != 7)
+ mexErrMsgTxt("expected (x1,x2,y1,y2,d1,d2,d3)");
+ if (nlhs > 2)
+ mexErrMsgTxt("has maximal two output arguments");
+
+
+ plhs[0] = mxCreateDoubleMatrix(1, 3, mxREAL);
+
+ *mxGetPr(plhs[0]) = calcParIntA(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+ *(mxGetPr(plhs[0])+1) = calcParIntQX1X2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+ *(mxGetPr(plhs[0])+2) = calcParIntQY1X1(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+
+}
+
+
+
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
- //printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-
- *mxGetPr(plhs[0]) = g0(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
- *mxGetPr(plhs[1]) = *mxGetPr(prhs[1])!=0?singleQuad(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3])):0;
-
+ *mxGetPr(plhs[0]) = g_AY(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+ *mxGetPr(plhs[0]) = g_QY(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
}
mex mex_Fpar.cpp slpRectangle.cpp
mex mex_Fort.cpp slpRectangle.cpp
-mex mex_FparY1Y2_2.cpp slpRectangle.cpp
-
mex mex_build_A.cpp slpRectangle.cpp
mex mex_build_As1.cpp slpRectangle.cpp
mex mex_build_As2.cpp slpRectangle.cpp
#include <cassert>\r
#include <stdlib.h>\r
\r
-#include <mex.h>\r
+//#include <mex.h>\r
\r
#include "slpRectangle.hpp"\r
\r
double gauss_w5[] = {0.1185,0.2393,0.2844,0.2393,0.1185};\r
double gauss_n5[] = {0.0469,0.2308,0.5000,0.7692,0.9531};\r
\r
+\r
+\r
int sign(double);\r
//double arsinh(double);\r
\r
+\r
+/* ============================== */\r
+\r
int inline sign(double x) {\r
return x > 0 ? 1 : (x < 0 ? -1 : 0);\r
}\r
}\r
*/\r
\r
-double singleQuad(double p, double y, double x, double l) {\r
+double inline max(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
\r
return sol;\r
}\r
\r
-\r
-double doubleQuad(double p, double y1, double y2, double x1, double x2, double l) {\r
-\r
- double sol = 0;\r
- for(int i = 0;i<5;++i){\r
- for(int j=0;j<5;++j){\r
- sol += pow((x1-y1*gauss_n5[i])*(x1-y1*gauss_n5[i])\r
- +(x2-y2*gauss_n5[j])*(x2-y2*gauss_n5[j])+l*l,p)*y1*gauss_w5[i]*y2*gauss_w5[j];\r
- }\r
- }\r
- return sol;\r
-}\r
-\r
//y-x muss != 0 sein\r
-double g0(double p, double y, double x, double l) {\r
+double g_AY(double p, double y, double x, double l) {\r
//printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);\r
\r
double sol = 0;\r
sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));\r
else\r
sol = (y - x) * pow((y - x) * (y - x) + l * l, p) + 2 * p * l * l\r
- * g0(p - 1, y, x, l) / (2 * p + 1);\r
+ * g_AY(p - 1, y, x, l) / (2 * p + 1);\r
} else {\r
if (p == -0.5)\r
sol = sign(y - x) * log(fabs(y - x));\r
return sol;\r
}\r
\r
-double G00(double p, double y1, double y2, double x1, double x2, double l) {\r
+double G_QY1Y2(double p, double y1, double y2, double x1, double x2, double l) {\r
+\r
+ double sol = 0;\r
+ for(int i = 0;i<5;++i){\r
+ for(int j=0;j<5;++j){\r
+ sol += pow((x1-y1*gauss_n5[i])*(x1-y1*gauss_n5[i])\r
+ +(x2-y2*gauss_n5[j])*(x2-y2*gauss_n5[j])+l*l,p)*y1*gauss_w5[i]*y2*gauss_w5[j];\r
+ }\r
+ }\r
+ return sol;\r
+}\r
+\r
+\r
+double G_AY2X2(double y1, double y2, double x1, double x2, double d3) {\r
+\r
+ // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
+\r
+ double hlp = ((x1 - y1 ) * (x1 - y1 ) + (x2 - y2) * (x2 - y2) + d3 * d3);\r
+\r
+ double sol = sqrt(hlp);\r
+\r
+ if ((x2 - y2) != 0)\r
+ sol += (x2 - y2) * log(sqrt(hlp)-(x2 - y2));\r
+\r
+ return sol;\r
+}\r
+\r
+double G_AY1Y2(double p, double y1, double y2, double x1, double x2, double l) {\r
// printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);\r
double pt = p + 1.5;\r
double sol = 0;\r
\r
} else if ((pt > 0) && ((int) pt == pt)) {\r
if (l != 0)\r
- sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l);\r
+ sol = 2 * p * l * l * G_AY1Y2(p - 1, y1, y2, x1, x2, l);\r
if ((y1 - x1) != 0)\r
- sol += (y1 - x1) * g0(p, y2, x2,\r
+ sol += (y1 - x1) * g_AY(p, y2, x2,\r
sqrt((y1 - x1) * (y1 - x1) + l * l));\r
if ((y2 - x2) != 0)\r
- sol += (y2 - x2) * g0(p, y1, x1,\r
+ sol += (y2 - x2) * g_AY(p, y1, x1,\r
sqrt((y2 - x2) * (y2 - x2) + l * l));\r
sol /= 2 * p + 2;\r
} else {\r
sol = NAN;\r
- cout << "warning in G00: no case for p="<< p <<" defined. Possible: [-1.5,-0.5,0.5,...]\n";\r
- //mexErrMsgTxt("no case for p defined\n");\r
+ cout << "warning in G_AY1Y2: no case for p="<< p <<" defined. Possible: [-1.5,-0.5,0.5,...]\n";\r
}\r
\r
return sol;\r
}\r
\r
-double F_parY1Y2(double x1, double x2, double y1, double y2, double d1, double d2,\r
- double d3) {\r
-\r
- // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
- double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
-\r
- if (sol != 0)\r
- sol *= G00(-0.5, x1, x2, y1 + d1, y2 + d2, d3);\r
-\r
- if ((x1 - y1 - d1) != 0)\r
- sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,\r
- sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));\r
-\r
- if ((x2 - y2 - d2) != 0)\r
- sol -= (x2 - y2 - d2) * g0(0.5, x2, y2 + d2,\r
- sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));\r
-\r
- double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
- - d2) + d3 * d3);\r
- sol += 1. / 3 * hlp * sqrt(hlp);\r
- return sol;\r
-}\r
-\r
-double F_parY1X2(double x1, double x2, double y1, double y2, double d1, double d2,\r
- double d3) {\r
-\r
- // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
- double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
-\r
- if (sol != 0)\r
- sol *= G00(-0.5, x1, -y2 - d2, y1 + d1, -x2, d3);\r
-\r
- if ((x1 - y1 - d1) != 0)\r
- sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,\r
- sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));\r
+double Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2, double l) {\r
+ double sol = 0;\r
\r
- if ((x2 - y2 - d2) != 0)\r
- sol -= (x2 - y2 - d2) * g0(0.5, -y2 - d2, -x2,\r
- sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));\r
+ //sol += y2*G_AY1Y2(y1,x2,)\r
\r
- double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
- - d2) + d3 * d3);\r
- sol += 1. / 3 * hlp * sqrt(hlp);\r
return sol;\r
}\r
\r
-double F_parX1Y2(double x1, double x2, double y1, double y2, double d1, double d2,\r
+double F_par(double x1, double x2, double y1, double y2, double d1, double d2,\r
double d3) {\r
\r
// printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
\r
if (sol != 0)\r
- sol *= G00(-0.5, -y1 - d1, x2, -x1, y2 + d2, d3);\r
+ sol *= G_AY1Y2(-0.5, x1, x2, y1 + d1, y2 + d2, d3);\r
\r
if ((x1 - y1 - d1) != 0)\r
- sol -= (x1 - y1 - d1) * g0(0.5, -y1 - d1, -x1,\r
+ sol -= (x1 - y1 - d1) * g_AY(0.5, x1, y1 + d1,\r
sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));\r
\r
if ((x2 - y2 - d2) != 0)\r
- sol -= (x2 - y2 - d2) * g0(0.5, x2, y2 + d2,\r
+ sol -= (x2 - y2 - d2) * g_AY(0.5, x2, y2 + d2,\r
sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));\r
\r
double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
return sol;\r
}\r
\r
-double F_parY1Y2_2(double x1, double x2, double y1, double y2, double d1, double d2,\r
- double d3) {\r
-\r
- // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
- double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
-\r
- if (sol != 0)\r
- sol *= doubleQuad(-0.5, x1, x2, y1 + d1, y2 + d2, d3);\r
-\r
- if ((x1 - y1 - d1) != 0)\r
- sol -= (x1 - y1 - d1) * singleQuad(0.5, x1, y1 + d1,\r
- sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));\r
\r
- if ((x2 - y2 - d2) != 0)\r
- sol -= (x2 - y2 - d2) * singleQuad(0.5, x2, y2 + d2,\r
- sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));\r
-\r
- double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
- - d2) + d3 * d3);\r
- sol += 1. / 3 * hlp * sqrt(hlp);\r
- return sol;\r
-}\r
-\r
-double F_parQY2X2(double x1, double x2, double y1, double y2, double d1, double d2,\r
- double d3) {\r
-\r
- // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
-\r
- double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
- - d2) + d3 * d3);\r
-\r
- double sol = sqrt(hlp);\r
-\r
- if ((x2 - y2 - d2) != 0)\r
- sol += (x2 - y2 - d2) * log(sqrt(hlp)-(x2 - y2 - d2));\r
-\r
- return sol;\r
-}\r
\r
double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,\r
double d3) {\r
\r
// printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y2,y3,d1,d2,d3);\r
- double sol = -G00(0.5, y3, x1, -d3, d1, x2 - y2 - d2);\r
+ double sol = -G_AY1Y2(0.5, y3, x1, -d3, d1, x2 - y2 - d2);\r
\r
if ((x1 - d1) * (x2 - y2 - d2) != 0)\r
- sol -= (x1 - d1) * (x2 - y2 - d2) * G00(-0.5, x2, y3, y2 + d2, -d3,\r
+ sol -= (x1 - d1) * (x2 - y2 - d2) * G_AY1Y2(-0.5, x2, y3, y2 + d2, -d3,\r
x1 - d1);\r
\r
if ((x1 - d1) != 0)\r
- sol += (x1 - d1) * g0(0.5, y3, -d3,\r
+ sol += (x1 - d1) * g_AY(0.5, y3, -d3,\r
sqrt((x1 - d1) * (x1 - d1) + (x2 - y2 - d2) * (x2 - y2 - d2)));\r
\r
if ((y3 + d3) * (x2 - y2 - d2) != 0)\r
- sol -= (y3 + d3) * (x2 - y2 - d2) * G00(-0.5, x1, x2, d1, y2 + d2,\r
+ sol -= (y3 + d3) * (x2 - y2 - d2) * G_AY1Y2(-0.5, x1, x2, d1, y2 + d2,\r
-y3 - d3);\r
\r
if ((y3 + d3) != 0)\r
- sol += (y3 + d3) * g0(0.5, x1, d1,\r
+ sol += (y3 + d3) * g_AY(0.5, x1, d1,\r
sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + (y3 + d3) * (y3 + d3)));\r
\r
return sol / 2.;\r
}\r
\r
-double applyInt4(\r
+/*double applyInt4(\r
double(*f)(double, double, double, double, double, double, double),\r
double a, double b, double c, double d, double r, double t, double u,\r
double v, double d1, double d2, double d3) {\r
d1, d2, d3) + f(a, d, r, v, d1, d2, d3) - f(a, d, r, u, d1, d2, d3)\r
+ f(a, c, t, v, d1, d2, d3) - f(a, c, t, u, d1, d2, d3) - f(a, c,\r
r, v, d1, d2, d3) + f(a, c, r, u, d1, d2, d3);\r
- /*\r
- return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(\r
- k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,\r
- s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,\r
- t1, l2, d1, d2, d3) - f(k1, s2, t1, t2, d1, d2, d3) - f(s1, k2, l1,\r
- l2, d1, d2, d3) + f(s1, k2, l1, t2, d1, d2, d3) + f(s1, k2, t1, l2,\r
- d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,\r
- d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,\r
- d3) + f(s1, s2, t1, t2, d1, d2, d3);\r
-\r
- */\r
-}\r
+\r
+}*/\r
\r
double apply0Int4(\r
double(*f)(double, double, double, double, double, double, double),\r
\r
}\r
\r
-double slpADLO(double y1, double y2, double x1, double x2, double a) {\r
- double G3 = 0;\r
- double gL = 0;\r
- double gK = 0;\r
- double tmp;\r
+// Berechnet das eigentliche Integral für parallele Elemente voll Analytisch\r
+double calcParIntA(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+ return apply0Int4(F_par, b, d, t, v, d1, d2, d3);\r
\r
- tmp = y1 - x1;\r
- if (fabs(tmp) >= EPS * y1) {\r
- tmp = sqrt(y1 * y1 + x1 * x1 + a * a - 2 * y1 * x1);\r
- gL = compute_g0(-0.5, y2, x2, tmp);\r
+}\r
+\r
+double calcParIntQX1X2(double b, double d, double t, double v, double d1,\r
+ double d2, double d3) {\r
+ //Fall 2\r
+ double sol = 0;\r
+\r
+ for (int i = 0; i < 5; ++i) {\r
+ for (int j = 0; j < 5; ++j) {\r
+ sol += G_AY1Y2(-0.5, t + d1, v + d2, b * gauss_n5[i],\r
+ d * gauss_n5[j], d3) * b * d * gauss_w5[i] * gauss_w5[j];\r
+ sol -= G_AY1Y2(-0.5, d1, v + d2, b * gauss_n5[i], d * gauss_n5[j],\r
+ d3) * b * d * gauss_w5[i] * gauss_w5[j];\r
+ sol -= G_AY1Y2(-0.5, t + d1, d2, b * gauss_n5[i], d * gauss_n5[j],\r
+ d3) * b * d * gauss_w5[i] * gauss_w5[j];\r
+ sol += G_AY1Y2(-0.5, d1, d2, b * gauss_n5[i], d * gauss_n5[j], d3)\r
+ * b * d * gauss_w5[i] * gauss_w5[j];\r
+ }\r
}\r
- tmp = y2 - x2;\r
- if (fabs(tmp) >= EPS * y2) {\r
- tmp = sqrt(y2 * y2 + x2 * x2 + a * a - 2 * y2 * x2);\r
- gK = compute_g0(-0.5, y1, x1, tmp);\r
+\r
+ return sol;\r
+\r
+}\r
+\r
+double calcParIntQY1X1(double b, double d, double t, double v, double d1,\r
+ double d2, double d3) {\r
+ //Fall 3\r
+ double sol = 0;\r
+\r
+ for (int i = 0; i < 5; ++i) {\r
+ for (int j = 0; j < 5; ++j) {\r
+ sol += G_AY2X2( t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i],d , d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
+ sol -= G_AY2X2( t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], 0, d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
+ sol -= G_AY2X2( t * gauss_n5[j] + d1, d2, b * gauss_n5[i], d, d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
+ sol += G_AY2X2( t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
+ }\r
}\r
- if (fabs(a * a) > EPS) {\r
- if ((y1 - x1) * (y2 - x2) * a >= 0)\r
- tmp = 1.;\r
- else\r
- tmp = -1.;\r
\r
- G3 = tmp * acos(\r
- (-2. * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2)) / (((y1\r
- - x1) * (y1 - x1) + a * a) * ((y2 - x2) * (y2 - x2) + a\r
- * a)) + 1.) / (2. * a);\r
+ return sol;\r
+\r
+}\r
+\r
+double calcParIntQY1(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+ //Fall 4\r
+ return apply0Int4(F_par, b, d, t, v, d1, d2, d3); ///ACHTUNG noch Falsch\r
+\r
+}\r
+\r
+double calcOrtIntA(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+ return apply0Int4(F_ort, b, d, t, v, d1, d2, d3);\r
+\r
+}\r
+\r
+double calcOrtIntQY1Y2(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+ return 0;\r
+\r
+}\r
+\r
+\r
+double calcParIntO1(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+\r
+ //Elmente anordnen\r
+ if(b*b+d*d>t*t+v*v){\r
+ double tmp = 0;\r
+\r
+ tmp = b; b = t; t = tmp;\r
+ tmp = d; d = v; v = tmp;\r
+\r
+ d1 = -d1;\r
+ d2 = -d2;\r
+ d3 = -d3;\r
}\r
\r
- return (y1 - x1) * gL + (y2 - x2) * gK - a * a * G3;\r
+ return calcParIntA( b, d, t, v, d1, d2, d3);\r
}\r
\r
-double compute_g0(double p, double y, double x, double a) {\r
- int sp = (int) 2 * (p - EPS); // MK (p-EPS) instead of (p+EPS)\r
- //printf("\n compute_g0, p = %lf, sp = %d\n",p,sp);\r
- assert(\r
- p == 0 || (p == -0.5) || (p == -1) || (p == -1.5) || (fabs(a)\r
- <= EPS));\r
- if (fabs(a) <= EPS) {\r
- // printf("\n a < eps\n");\r
- switch (sp) {\r
- case 0:\r
- return y - x;\r
- case -1:\r
- return log(fabs(x - y)) * (y - x) / fabs(y - x);\r
- case -2:\r
- return -(y - x) / fabs(y - x) / fabs(y - x);\r
- case -3:\r
- return -0.5 * (y - x) / fabs(y - x) / fabs(y - x) / fabs(y - x);\r
- }\r
- } else {\r
- // printf("\n a > eps\n");\r
- switch (sp) {\r
- case 0:\r
- return y - x;\r
- case -1:\r
- return asinh((y - x) / fabs(a));\r
- case -2:\r
- return atan((y - x) / fabs(a));\r
- case -3:\r
- return (y - x) * pow((x * x + y * y + a * a - 2 * x * y), -0.5)\r
- / (a * a);\r
- default:\r
- cout << "p must be either 0, -1/2, -1 or -3/2. (" << p << ")\n";\r
- return NAN;\r
- }\r
+double calcOrtIntO1(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+\r
+ //Elmente anordnen\r
+ if(b*b+d*d>t*t+v*v){\r
+ double tmp = 0;\r
+\r
+ tmp = b; b = t; t = tmp;\r
+ tmp = d; d = v; v = tmp;\r
+\r
+ d1 = -d1;\r
+ d2 = -d2;\r
+ d3 = -d3;\r
+ }\r
+\r
+ return calcOrtIntA( b, d, t, v, d1, d2, d3);\r
+}\r
+\r
+\r
+double calcParIntO2(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {\r
+\r
+ //Elmente anordnen\r
+ if(b*b+d*d>t*t+v*v){\r
+ double tmp = 0;\r
+\r
+ tmp = b; b = t; t = tmp;\r
+ tmp = d; d = v; v = tmp;\r
+\r
+ d1 = -d1;\r
+ d2 = -d2;\r
+ d3 = -d3;\r
+ }\r
+\r
+ if((b*b+d*d) *eta <= d1*d1+d2*d2+d3*d3){\r
+ return calcParIntQX1X2( b, d, t, v, d1, d2, d3);\r
+ }else{\r
+ return calcParIntA( b, d, t, v, d1, d2, d3);\r
+ }\r
+\r
+\r
+}\r
+\r
+double calcOrtIntO2(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {\r
+\r
+ //Elmente anordnen\r
+ if(b*b+d*d>t*t+v*v){\r
+ double tmp = 0;\r
+\r
+ tmp = b; b = t; t = tmp;\r
+ tmp = d; d = v; v = tmp;\r
+\r
+ d1 = -d1;\r
+ d2 = -d2;\r
+ d3 = -d3;\r
+ }\r
+\r
+ if((b*b+d*d) *eta <= d1*d1+d2*d2+d3*d3){\r
+ return 0;\r
+ }else{\r
+ return calcOrtIntA( b, d, t, v, d1, d2, d3);\r
}\r
}\r
\r
-double FLO_plane(double x1, double x2, double y1, double y2, double delta1,\r
- double delta2, double a) {\r
- double yd1 = y1 + delta1;\r
- double yd2 = y2 + delta2;\r
- double tmp1 = x1 - y1 - delta1;\r
- double tmp2 = x2 - y2 - delta2;\r
- double tmp3 = sqrt(tmp2 * tmp2 + a * a);\r
- double tmp4 = sqrt(tmp1 * tmp1 + a * a);\r
- double tmp5 = pow(tmp1 * tmp1 + tmp2 * tmp2 + a * a, 3. / 2.);\r
- double rval = 0;\r
-\r
- rval = tmp1 * tmp2 * slpADLO(x1, x2, yd1, yd2, a) - tmp1 * compute_g0(-0.5,\r
- x1, yd1, tmp3) - tmp2 * compute_g0(-0.5, x2, yd2, tmp4) + tmp5 / 3.;\r
- return rval;\r
+\r
+double calcParIntO3(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {\r
+\r
+ //Achsen vertauschen\r
+ if(b*b+t*t>d*d+v*v){\r
+ double tmp = 0;\r
+\r
+ tmp = b; b = d; d = tmp;\r
+ tmp = t; t = v; v = tmp;\r
+\r
+ tmp = d1; d1 = d2; d2 = tmp;\r
+ }\r
+\r
+ if(max(b,t) *eta <= d1){\r
+ return calcParIntQY1X1(b, d, t, v, d1, d2, d3);\r
+ }else{\r
+ return calcParIntA( b, d, t, v, d1, d2, d3);\r
+ }\r
+\r
+\r
}\r
\r
-double FLO_perp(double x1, double x2, double y2, double y3, double delta1,\r
- double delta2, double a) {\r
- double xd1 = x1 - delta1;\r
- double xd2 = x2 - delta2;\r
- double yd2 = y2 - delta2;\r
- double yd3 = y3 - a;\r
- double tmp2 = x2 - y2 - delta2;\r
- double rval = 0;\r
-\r
- rval = xd1 * compute_g0(y3, a, sqrt(xd1 * xd1 + tmp2 * tmp2), -0.5) + yd3\r
- * compute_g0(x1, delta1, sqrt(tmp2 * tmp2 + yd3 * yd3), -0.5) - xd1\r
- * tmp2 * slpADLO(x2, y3, yd2, a, xd1) - yd3 * tmp2 * slpADLO(x1,\r
- x2, delta1, yd2, -y3 - a) - slpADLO(y3, x1, -a, delta1, tmp2);\r
- return rval;\r
+double calcOrtIntO3(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {\r
+\r
+ //Elmente anordnen\r
+ if(b*b+d*d>t*t+v*v){\r
+ double tmp = 0;\r
+\r
+ tmp = b; b = t; t = tmp;\r
+ tmp = d; d = v; v = 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
+ return calcOrtIntA( b, d, t, v, d1, d2, d3);\r
+ }\r
}\r
#ifndef HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
#define HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
-int sign(double);
-double arsinh(double);
+//int sign(double);
+//double arsinh(double);
-// sol = g0(p,y,x,l); Quadratur !!! 0-y
-double singleQuad(double, double, double, double);
-
-// sol = G00(p,y1,y2,x1,x2,l); Quadratur !!! 0-y1 & 0-y2
-double doubleQuad(double, double, double, double, double, double);
// sol = g0(p,y,x,l);
-double g0(double, double, double, double);
-// sol = G00(p,y1,y2,x1,x2,l);
-double G00(double, double, double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-//double F_par(double, double, double, double, double, double, double);
-#define F_par F_parY1Y2
-// sol = F_ort(x1,x2,y2,y3,d1,d2,d3);
-double F_ort(double, double, double, double, double, double, double);
-
-
+double g_AY(double, double, double, double);
+// sol = g0(p,y,x,l);
+double g_QY(double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-double F_parY1Y2(double, double, double, double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-double F_parY1X2(double, double, double, double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-double F_parX1Y2(double, double, double, double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3); Quadratur!!!
-double F_parY1Y2_2(double, double, double, double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3); Quadratur!!!
-double F_parQY2X2(double, double, double, double, double, double, double);
+// sol = G00(p,y1,y2,x1,x2,l);
+double G_AY1Y2(double, double, double, double, double, double);
+// sol = G00(p,y1,y2,x1,x2,l);
+double G_AY2X2(double, double, double, double, double, double);
+// sol = G00(p,y1,y2,x1,x2,l);
+double G_QY1Y2(double, double, double, double, double, double);
-// funktionen von Maik
-double slpADLO(double, double, double, double, double);
-double compute_g0(double, double, double, double);
-double FLO_plane(double, double, double, double, double, double, double);
-double FLO_perp(double, double, double, double, double, double, double);
+// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
+double F_par(double, double, double, double, double, double, double);
+// sol = F_ort(x1,x2,y2,y3,d1,d2,d3);
+double F_ort(double, double, double, double, double, double, double);
-// sol = quadInt((F_par/F_ort),a,b,c,d,r,t,u,v,d1,d2,d3);
-double applyInt4(
- double(*)(double, double, double, double, double, double, double),
- double, double, double, double, double, double, double, double, double,
- double, double);
-// sol = quad0Int((F_par/F_ort),b,d,t,v,d1,d2,d3);
+// sol = quad0Int4((F_par/F_ort),b,d,t,v,d1,d2,d3);
double apply0Int4(
double(*)(double, double, double, double, double, double, double),
double, double, double, double, double, double, double);
-double apply0Int2(
- double(*)(double, double, double, double, double, double, double),
- double, double, double, double, double, double, double);
+// sol = calcParInt.(b,d,t,v,d1,d2,d3);
+double calcParIntA(double, double, double, double, double, double, double);
+double calcParIntQX1X2(double, double, double, double, double, double, double);
+double calcParIntQY1X1(double, double, double, double, double, double, double);
+double calcParIntQX1(double, double, double, double, double, double, double);
-double calcParIntA(double, double, double, double, double, double, double);
+// sol = calcOrtInt.(b,d,t,v,d1,d2,d3);
+double calcOrtIntA(double, double, double, double, double, double, double);
+double calcOrtIntQX1X2(double, double, double, double, double, double, double);
+// Elemente Vertauschen
+// sol = calc...Int01(b,d,t,v,d1,d2,d3);
+double calcParIntO1(double, double, double, double, double, double, double);
+double calcOrtIntO1(double, double, double, double, double, double, double);
+// Quadratur ueber kleineres Element
+// sol = calc...Int02(b,d,t,v,d1,d2,d3,eta);
+double calcParIntO2(double, double, double, double, double, double, double, double);
+double calcOrtIntO2(double, double, double, double, double, double, double, double);
+// Quadratur kuerzere Achse/Seite
+// sol = calc...Int03(b,d,t,v,d1,d2,d3,eta);
+double calcParIntO3(double, double, double, double, double, double, double, double);
+double calcOrtIntO3(double, double, double, double, double, double, double, double);
#endif
+++ /dev/null
-/*
- * test.cpp
- *
- * Created on: 06.04.2011
- * Author: treecity
- */
-
-#include <iostream>
-#include <cmath>
-
-using namespace std;
-
-
- int static sign(double x)
- {
- return x > 0 ? 1 : (x < 0 ? -1 : 0);
- }
-
- /***
- * ! Diese Funktionen benutzen nur die ersten beiden Dimensionen !
- */
-
- double g(double p, double y, double x, double lambda){
- double sol = NAN;
-
- if(lambda==0){ //Lemma B1
- if(p!=-0.5)
- sol = (y-x)*pow(fabs(y-x),2*p)/(2*p+1); //(p<.5) & x außerhalb ?
- else
- sol = sign(y-x)*log(fabs(y-x)); //x außerhalb ?
- }else{ //Lemma B2 (wird eigententlich nicht gebraucht)
-
-
- }
- return sol;
- }
-
- double G(double p, double * y, double * x, double lambda) {
- double sol =NAN;
- if(p==-1.5){
- if(lambda==0){
-
- }else{
-
- }
-
- }else if(p==-0.5){
-
- }
-
-
- return sol;
- }
-
- double F(double* x, double * y, double delta1, double delta2,double delta3){
-
- return NAN;
- }
-
-int main(){
-
- double p=-2,y=2,x=3.5,l=0;
-
-
- cout << "g(" << p << "," << y << "," << x << "," << l << ")=" << g(p,y,x,l) << "\n";
- p = -0.5;
- cout << "g(" << p << "," << y << "," << x << "," << l << ")=" << g(p,y,x,l) << "\n";
- x = 1.5;
- cout << "g(" << p << "," << y << "," << x << "," << l << ")=" << g(p,y,x,l) << "\n";
- p = 0;
- cout << "g(" << p << "," << y << "," << x << "," << l << ")=" << g(p,y,x,l) << "\n";
- l = -2;
- cout << "g(" << p << "," << y << "," << x << "," << l << ")=" << g(p,y,x,l) << "\n";
-
-
-
-
- return 0;
-}
sol3 = surfQuad(@(x1,y1)(int2(x1,y1,c,d,u,v)),a,b,r,t,5);
- ret = [quad_sol sol sol2 sol3];
+ sol4 =...
+ surfQuad(@(x1,y1)mex_FparQY2X2(x1,d,y1,v,l1,l2,l3),a,b,r,t,5)...
+ -surfQuad(@(x1,y1)mex_FparQY2X2(x1,c,y1,v,l1,l2,l3),a,b,r,t,5)...
+ -surfQuad(@(x1,y1)mex_FparQY2X2(x1,d,y1,u,l1,l2,l3),a,b,r,t,5)...
+ +surfQuad(@(x1,y1)mex_FparQY2X2(x1,c,y1,u,l1,l2,l3),a,b,r,t,5);
+
+
+ ret = [quad_sol sol sol2 sol3 sol4];
end
--- /dev/null
+\r
+mex mex_build_As3.cpp slpRectangle.cpp\r
+\r
+dataIso =[];\r
+dataAniso =[];\r
+maxSize = 500;\r
+fileName = './test_save';\r
+Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren\r
+% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren\r
+% Netz = 'exmpl_2DQuad';\r
+\r
+\r
+%% Anisotrope Verfeinerung\r
+\r
+if(~exist('d1'))\r
+ \r
+ disp 'Normal'\r
+ load(Netz)\r
+\r
+ ref = 0;\r
+\r
+while size(elements,1)<maxSize\r
+ tic\r
+ ref = ref+1; \r
+ \r
+\r
+\r
+ %Netz Isotrop verfeinern (4 Teile)\r
+ [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
+ \r
+ %Matrix mit feinem Netz berechnen\r
+ A_fine = mex_build_A(coordinates_fine,elements_fine);\r
+ \r
+ %rechte Seite Aufbauen\r
+ b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
+ \r
+ %Lösungsvektor bestimmen\r
+ x_fine = A_fine\b;\r
+ \r
+ %Energienorm^2 bestimmen\r
+ xe_fine = x_fine'*A_fine*x_fine;\r
+\r
+ % Welche Elemente sollen verfeinert werden\r
+ ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
+ \r
+ % Wie soll verfeinert werden\r
+ marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
+\r
+ %Daten zur Auswertung zwischenspeichern (testweise?)\r
+ dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
+ \r
+ %Sicherheit falls keine Elemente markiert sind\r
+ if(sum(marked~=1)==0)\r
+ disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
+ break;\r
+ end\r
+ \r
+ %Verfeinern\r
+ %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
+ file = ['save' int2str(ref)];\r
+ load(file);\r
+ toc\r
+end\r
+d1 = dataAniso\r
+end\r
+\r
+%% Anisotrope Verfeinerung\r
+disp 'Optimiert'\r
+load(Netz)\r
+\r
+ref = 0;\r
+while size(elements,1)<maxSize\r
+ tic\r
+ ref = ref+1; \r
+ \r
+\r
+\r
+ %Netz Isotrop verfeinern (4 Teile)\r
+ [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
+ \r
+ %Matrix mit feinem Netz berechnen\r
+ A_fine = mex_build_As3(coordinates_fine,elements_fine);\r
+ \r
+ %rechte Seite Aufbauen\r
+ b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
+ \r
+ %Lösungsvektor bestimmen\r
+ x_fine = A_fine\b;\r
+ \r
+ %Energienorm^2 bestimmen\r
+ xe_fine = x_fine'*A_fine*x_fine;\r
+\r
+ % Welche Elemente sollen verfeinert werden\r
+ ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
+ \r
+ % Wie soll verfeinert werden\r
+ marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
+\r
+ %Daten zur Auswertung zwischenspeichern (testweise?)\r
+ dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
+ \r
+ %Sicherheit falls keine Elemente markiert sind\r
+% if(sum(marked~=1)==0)\r
+% disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
+% break;\r
+% end\r
+ \r
+ %Verfeinern\r
+ %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
+ file = ['save' int2str(ref)];\r
+ load(file);\r
+ toc\r
+end\r
+d2 = dataAniso\r
+\r
+clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine\r
+\r
+\r
+%% Alles Zeichnen\r
+figure (1)\r
+loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));\r
+\r
+legend('Normal','optimiert')\r
+xlabel('Elemente')\r
+ylabel('Fehler')\r
+\r
+figure (2)\r
+loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));\r
+\r
+legend('Normal','Optimiert')\r
+xlabel('Elemente')\r
+ylabel('EnergieNorm')\r
projects / bacc.git / commitdiff