leg1 = {};
leg2 = {};
leg3 = {};
+leg4 = {};
+
for i = 1:length(files)
load(files{i});
[m n] = size(data);
- step = round(n/3);
+ step = round(n/4);
if step<1
disp (['Error: No Data to show. : ' i])
continue;
end
G_D = [G_D data(:,2:end)];
- X = [X repmat(data(:,1),1,step*3)];
+ X = [X repmat(data(:,1),1,step*4)];
p1 = find(files{i}=='/',1,'last')+1;
p2 = find(files{i}(p1:end)=='_',1);
l0 = [files{i}(p1:p1+p2-2) ' '];
- l1 = {type2str{data(1,[2+(0:step-1)*3])}};
+ l1 = {type2str{data(1,[2+(0:step-1)*4])}};
l2 ={};
l3 ={};
l4 ={};
for i = 1:step
l2{i} = [ l0 l1{i}];
- l3{i} = ['\mu ' l0 l1{i}];
- l4{i} = ['error ' l0 l1{i}];
+ l3{i} = ['C^{-1}\mu ' l0 l1{i}];
+ l4{i} = ['\mu ' l0 l1{i}];
+ l5{i} = ['error ' l0 l1{i}];
end
leg1 = {leg1{:} l2{:}}';
leg2 = {leg2{:} l3{:}}';
leg3 = {leg3{:} l4{:}}';
+ leg4 = {leg4{:} l5{:}}';
end
disp ('Error: No Data to show.')
else
- data((end-8):end,[1 3 4])
+ data((end-8):end,[1 3 4 5])
sol = interp1(1./X((1):(end-2),3)',G_D((1):(end-2),3)',0,'spline')
% sol = 8.28466;
figure(4)
loglog(X(:,[2+(0:step-1)*3]),G_D(:,[2+(0:step-1)*3]),'--o')
hold on
+loglog(X(:,[3+(0:step-1)*3]), G_D(:,3+(0:step-1)*3)*G_D(1,2)/G_D(1,3),':x')
+loglog(X(:,[4+(0:step-1)*3]), sqrt(abs(sol -G_D(:,4+(0:step-1)*3)))*G_D(1,2)/sqrt(abs(sol-G_D(1,4))),'-*')
loglog(X(:,1),[20*X(:,1).^(-1/2),6*X(:,1).^(-1/4),5*X(:,1).^(-3/4)],'-.')
-loglog(X(:,[3+(0:step-1)*3]), sqrt(abs(sol -G_D(:,3+(0:step-1)*3)))*G_D(1,2)/sqrt(abs(sol-G_D(1,3))),'-*')
+
hold off
title('Fehler')
xlabel('Elemente');
ylabel('Schaetzer');
-legend({leg2{:}...
+legend({leg2{:} leg3{:} leg4{:}...
'N^{-1/2}' 'N^{-1/4}' 'N^{-3/4}'...
- leg3{:}} ,'location','southwest','box', 'off');
+ } ,'location','southwest','box', 'off');
print('-r600','-depsc',[printt '_error.eps'])
figure(5)
-loglog(X(:,[3+(0:step-1)*3]),G_D(:,[3+(0:step-1)*3]),'-*',X(:,1),repmat(sol,size(X,1),1),'-.')
-ylim([min(min(G_D(:,[3+(0:step-1)*3]))) 1.005*max(max(G_D(:,[3+(0:step-1)*3])))])
+loglog(X(:,[4+(0:step-1)*3]),G_D(:,[4+(0:step-1)*3]),'-*',X(:,1),repmat(sol,size(X,1),1),'-.')
+ylim([min(min(G_D(:,[4+(0:step-1)*3]))) 1.005*max(max(G_D(:,[4+(0:step-1)*3])))])
title('Energie Norm')
xlabel('Elemente');
ylabel('eNorm^2');
%benötigte Zeit Schätzen
if(size(G_T,1)>2)
nextF = G_T(end,1)*4;
- nextS = neville2(1:size(G_T,1),G_T(:,1)',size(G_T,1)+1);
+ nextS = interp1(1:size(G_T,1),G_T(:,1)',size(G_T,1)+1,'slpine');
cald = size(G_T,1)+(-1:0);
calc = size(G_T,1)+(-2:0);
- nextTime = [nextS neville2(G_T(calc,1)',G_T(calc,2)',nextF) ...
- neville2(G_T(cald,1)',G_T(cald,3)',nextF) ...
- neville2(G_T(calc,1)',G_T(calc,4)',nextS)];
+ nextTime = [nextS interp1(G_T(calc,1)',G_T(calc,2)',nextF,'spline') ...
+ interp1(G_T(cald,1)',G_T(cald,3)',nextF,'spline') ...
+ interp1(G_T(calc,1)',G_T(calc,4)',nextS,'spline')];
end
% Ein kompletter Verfeinerungschritt
%Lösung Berechnen
x_fine = A_fine\b;
- %Energienorm^2 Berechnen
- xe_fine = x_fine'*A_fine*x_fine;
%Fehlerschätzer aufbauen
ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
+ %Fehlerschätzer 2 aufbauen
+ A = mex_build_AU(G_C,G_E,mu,type(i));
+ x = A\(sqrt(sum(quadNorm(G_C,G_E,'w').^2,2)));
+ xo_fine(f2s) = repmat(x,1,4);
+ xd_fine = xo_fine'-x_fine;
+ ind2 = xd_fine'*A_fine*xd_fine;
+
+
+ %Energienorm^2 Berechnen
+ xe_fine = x_fine'*A_fine*x_fine;
+% xe_fine = x'*A*x;
+
+% sqrt(sum(ind))
+% ind2
+
save_A{i} = A_fine;
save_x{i} = x_fine;
- data = [data type(i) sqrt(sum(ind)) xe_fine];
+
+
+
+ data = [data type(i) sqrt(sum(ind)) ind2 xe_fine];
end
time(2) = toc;
--- /dev/null
+function export_mesh(coo, ele, nei, f2s, file)
+
+plotShape(coo,ele,'tb');view(2);
+print('-r600','-depsc',['../doc/fig/' file '_ref.eps'])
+
+%% Koordinaten
+fid = fopen(['../doc/fig/' file '_coo.tex'],'w');
+% fprintf(fid,'\\begin{figure}\n');
+fprintf(fid,'\\begin{tabular}{>{\\columncolor{gray}}rccc}\n');
+fprintf(fid,'\\rowcolor{gray}\n Index & x1 & x2 & x3\\\\');
+
+[m n] =size(coo);
+
+for i = 1:m
+ if(i~=1)
+ fprintf(fid,'\\\\');
+ end
+ fprintf(fid,'\n %d & %.4f & %.4f & %.4f',i,coo(i,1),coo(i,2),coo(i,3));
+end
+fprintf(fid,'\n\\end{tabular}');
+% fprintf(fid,['\n\\caption{' file ' - Koordinaten}']);
+% fprintf(fid,['\n\\label{' file ':coo}']);
+% fprintf(fid,'\n\\end{figure}');
+fprintf(fid,'\n');
+fclose(fid);
+
+%% Elemente
+fid = fopen(['../doc/fig/' file '_ele.tex'],'w');
+% fprintf(fid,'\\begin{figure}\n');
+fprintf(fid,'\\begin{tabular}{>{\\columncolor{gray}}rcccc}\n');
+fprintf(fid,'\\rowcolor{gray}\n Index & c1 & c2 & c3 & c4\\\\');
+
+[m n] =size(ele);
+
+for i = 1:m
+ if(i~=1)
+ fprintf(fid,'\\\\');
+ end
+ fprintf(fid,'\n %d & %d & %d & %d & %d',i,ele(i,1),ele(i,2),ele(i,3),ele(i,4));
+end
+fprintf(fid,'\n\\end{tabular}');
+% fprintf(fid,['\n\\caption{' file ' - Elemente}']);
+fprintf(fid,['\n\\label{' file ':ele}']);
+% fprintf(fid,'\n\\end{figure}');
+fclose(fid);
+
+%% Nachbarn
+fid = fopen(['../doc/fig/' file '_nei.tex'],'w');
+% fprintf(fid,'\\begin{figure}\n');
+fprintf(fid,'\\begin{tabular}{>{\\columncolor{gray}}rcccccccc}\n');
+fprintf(fid,'\\rowcolor{gray}\n Index & n1 & n2 & n3 & n4 & n5 & n6 & n7 & n8\\\\');
+
+[m n] =size(nei);
+
+for i = 1:m
+ if(i~=1)
+ fprintf(fid,'\\\\');
+ end
+ fprintf(fid,'\n %d & %d & %d & %d & %d & %d & %d & %d & %d',i,nei(i,1),nei(i,2),nei(i,3),nei(i,4),nei(i,5),nei(i,6),nei(i,7),nei(i,8));
+end
+fprintf(fid,'\n\\end{tabular}');
+% fprintf(fid,['\n\\caption{' file ' - Nachbarn}']);
+fprintf(fid,['\n\\label{' file ':nei}']);
+% fprintf(fid,'\n\\end{figure}');
+fclose(fid);
+
+%% VaterSohn
+fid = fopen(['../doc/fig/' file '_f2s.tex'],'w');
+% fprintf(fid,'\\begin{figure}\n');
+fprintf(fid,'\\begin{tabular}{>{\\columncolor{gray}}rcccc}\n');
+fprintf(fid,'\\rowcolor{gray}\n Index & e1 & e2 & e3 & e4\\\\');
+
+[m n] =size(f2s);
+
+for i = 1:m
+ if(i~=1)
+ fprintf(fid,'\\\\');
+ end
+ fprintf(fid,'\n %d & %d & %d & %d & %d',i,f2s(i,1),f2s(i,2),f2s(i,3),f2s(i,4));
+end
+fprintf(fid,'\n\\end{tabular}');
+% fprintf(fid,['\n\\caption{' file ' - VaterSohn}']);
+fprintf(fid,['\n\\label{' file ':f2s}']);
+% fprintf(fid,'\n\\end{figure}');
+fclose(fid);
+end
+
*\r
*\r
*/\r
-\r
-#include <iostream>\r
+//#include <iostream>\r
#include <cmath>\r
-#include <cassert>\r
+//#include <cassert>\r
#include "mex.h"\r
-#include <stdlib.h>\r
+/*/#include <stdlib.h>\r
\r
//#include "gauss.hpp"\r
-\r
+*/\r
#define M_EPS 10e-8\r
-\r
+/*\r
//#include "tbb/parallel_for.h"\r
-\r
+*/\r
#include "slpRectangle.hpp"\r
-\r
-using namespace std;\r
-using namespace slpR;\r
-\r
+/*\r
+//using namespace std;\r
+//using namespace slpR;\r
+*/\r
int dimOfVec(double* vec) {\r
- std::cout << "";\r
if (vec[2] != 0)\r
return 2;\r
else if (vec[1] != 0)\r
return 1;\r
else if (vec[0] != 0)\r
return 0;\r
- else {\r
- cerr << "dimOfVec : (" << vec[0] << " " << vec[1] << " " << vec[2]\r
- << ") alle Werte 0 \n";\r
+ else {/*\r
+// cerr << "dimOfVec : (" << vec[0] << " " << vec[1] << " " << vec[2]\r
+// << ") alle Werte 0 \n";*/\r
return -1;\r
}\r
\r
}\r
\r
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
-\r
+/*\r
// initeQuad();\r
// cout << Quad[1].nod[0];\r
-\r
+*/\r
int i, j, k; //Schleifenindizes\r
double tmp; //Zwischenspeicherung der Einzelnen Werte\r
int count;\r
int cn = mxGetN(prhs[0]);\r
if (cn != 3)\r
mexErrMsgTxt("expected coordinates (Nx3)");\r
- cout << " Koordinaten:" << cm << endl;\r
+// cout << " Koordinaten:" << cm << endl;\r
\r
int em = mxGetM(prhs[1]);\r
int en = mxGetN(prhs[1]);\r
if (en != 4)\r
mexErrMsgTxt("expected elements (Mx4)");\r
- cout << " Elemente:" << em << endl;\r
+// cout << " Elemente:" << em << endl;\r
\r
//Auslesen der Parameter\r
\r
double*);\r
\r
//Art der Berechnung bestimmen\r
- cout << " Typ:" << type << endl;\r
+// cout << " Typ:" << type << endl;\r
switch (type) {\r
default:\r
ctypeP = cParO1;\r
}\r
\r
}\r
- cout << endl;\r
+// cout << endl;\r
//Rueckgabe (eventuell zurueck schreiben)\r
/* mxFree(x0);\r
mxFree(x1);\r
if(t)
for idx = eles
current = sum(coordinates(elements(idx,[2,4])',:),1)/2;
+ if(e)
+ cola = 'w';
+ else
+ cola = 'bla';
+ end
if(isempty(desc))
- text(current(1),current(2),current(3),num2str(idx),'color','w');
+ text(current(1),current(2),current(3),num2str(idx),'color',cola);
else
- text(current(1),current(2),current(3),desc{idx},'color','w');
+ text(current(1),current(2),current(3),desc{idx},'color',cola);
end
hold on
end
-#include <iostream>\r
-#include <cmath>\r
-#include <cassert>\r
-#include <limits>\r
-#include <cstdlib>\r
+//#include <iostream>\r
+#include <math.h>\r
+//#include <cassert>\r
+//#include <limits>\r
+//#include <cstdlib>\r
\r
//#include <mex.h>\r
\r
\r
#define quad 4 // Anzahl der Quadratur Auswertungstellen 2^quad\r
\r
-using namespace std;\r
+//using namespace std;\r
\r
int sign(double);\r
//double arsinh(double);\r
//sol -= (x2 - y2) * g_AY(-0.5,x2,y2, sqrt((x1 - y1) * (x1 - y1) + d3 * d3));\r
sol += (x2 - y2) * log(sqrt(hlp) - (x2 - y2));\r
\r
- if (sol != sol || fabs(sol) == numeric_limits<double>::infinity()) {\r
- cout << "G_AY2X2: " << sol << " isn't a Number. (" << y1 << "," << y2\r
- << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl;\r
- cout << (x2 - y2) << "__" << sqrt(hlp) << endl;\r
- cout << (x1 - y1) * (x1 - y1) << " " << hlp << endl;\r
+ if (sol != sol /*|| fabs(sol) == numeric_limits<double>::infinity()*/) {\r
+ // cout << "G_AY2X2: " << sol << " isn't a Number. (" << y1 << "," << y2\r
+ // << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl;\r
+ // cout << (x2 - y2) << "__" << sqrt(hlp) << endl;\r
+ // cout << (x1 - y1) * (x1 - y1) << " " << hlp << endl;\r
}\r
return sol;\r
}\r
sol /= 2 * p + 2;\r
} else {\r
sol = NAN;\r
- cout << "warning in G_AY1Y2: no case for p=" << p\r
- << " defined. Possible: [-1.5,-0.5,0.5,...]" << endl;\r
+ // cout << "warning in G_AY1Y2: no case for p=" << p\r
+ // << " defined. Possible: [-1.5,-0.5,0.5,...]" << endl;\r
}\r
\r
return sol;\r
}\r
\r
// Berechnet das eigentliche Integral fuer parallele Elemente voll Analytisch\r
-double slpR::calcParIntA(double b, double d, double t, double v, double d1,\r
+double calcParIntA(double b, double d, double t, double v, double d1,\r
double d2, double d3) {\r
return apply0Int4(F_par, b, d, t, v, d1, d2, d3);\r
\r
}\r
\r
-double slpR::calcParIntQX1X2(double b, double d, double t, double v, double d1,\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
}\r
\r
-double slpR::calcParIntQY1X1(double b, double d, double t, double v, double d1,\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
b * gauss_nodes[quad][i].n, 0, d3) * t * b\r
* gauss_nodes[quad][i].w * gauss_nodes[quad][j].w;\r
\r
- if (sol != sol || fabs(sol) == numeric_limits<double>::infinity()) {\r
- cout << "->(" << i << "," << j << ")" << endl;\r
+ if (sol != sol /*|| fabs(sol) == numeric_limits<double>::infinity()*/) {\r
+/* cout << "->(" << i << "," << j << ")" << endl;\r
cout << "-|("\r
<< G_AY2X2(t * gauss_nodes[quad][j].n + d1, v + d2,\r
b * gauss_nodes[quad][i].n, d, d3) * t * b\r
<< t * b * gauss_nodes[quad][i].w\r
* gauss_nodes[quad][j].w << ")" << endl;\r
\r
- cout << endl;\r
+ cout << endl;*/\r
return sol;\r
}\r
\r
}\r
}\r
\r
- if (sol != sol || sol == numeric_limits<double>::infinity())\r
- cout << "calcParIntQY1X1: " << sol << " isn't a Number. (" << b << ","\r
- << d << ")(" << t << "," << v << ")(" << d1 << "," << d2 << ","\r
- << d3 << ")" << endl;\r
+ //if (sol != sol /*|| sol == numeric_limits<double>::infinity()*/)\r
+ // cout << "calcParIntQY1X1: " << sol << " isn't a Number. (" << b << ","\r
+ // << d << ")(" << t << "," << v << ")(" << d1 << "," << d2 << ","\r
+ // << d3 << ")" << endl;\r
\r
return sol;\r
}\r
\r
-double slpR::calcParIntQY1(double b, double d, double t, double v, double d1,\r
+double calcParIntQY1(double b, double d, double t, double v, double d1,\r
double d2, double d3) {\r
//Fall 4\r
\r
\r
}\r
\r
-double slpR::calcParIntQ(double b, double d, double t, double v, double d1,\r
+double calcParIntQ(double b, double d, double t, double v, double d1,\r
double d2, double d3) {\r
//Fall 0\r
\r
return sol;\r
}\r
\r
-double slpR::calcOrtIntA(double b, double d, double t, double v, double d1,\r
+double calcOrtIntA(double b, double d, double t, double v, double d1,\r
double d2, double d3) {\r
return apply0Int4(F_ort, b, d, t, v, d1, d2, d3);\r
\r
}\r
\r
-double slpR::calcOrtIntQX1X2(double b, double d, double t, double v, double d1,\r
+double calcOrtIntQX1X2(double b, double d, double t, double v, double d1,\r
double d2, double d3) {\r
//Fall 2\r
double sol = 0;\r
\r
}\r
\r
-double slpR::calcParIntO0(double b, double d, double t, double v, double d1,\r
+double calcParIntO0(double b, double d, double t, double v, double d1,\r
double d2, double d3, double eta) {\r
return calcParIntA(b, d, t, v, d1, d2, d3);\r
}\r
-double slpR::calcOrtIntO0(double b, double d, double t, double v, double d1,\r
+double calcOrtIntO0(double b, double d, double t, double v, double d1,\r
double d2, double d3, double eta) {\r
return calcOrtIntA(b, d, t, v, d1, d2, d3);\r
}\r
\r
-double slpR::calcParIntO1(double b, double d, double t, double v, double d1,\r
+double calcParIntO1(double b, double d, double t, double v, double d1,\r
double d2, double d3, double eta) {\r
\r
//Elmente anordnen\r
return calcParIntA(b, d, t, v, d1, d2, d3);\r
}\r
\r
-double slpR::calcOrtIntO1(double b, double d, double t, double v, double d1,\r
+double calcOrtIntO1(double b, double d, double t, double v, double d1,\r
double d2, double d3, double eta) {\r
\r
//Elmente anordnen\r
return calcOrtIntA(b, d, t, v, d1, d2, d3);\r
}\r
\r
-double slpR::calcParIntO2(double b, double d, double t, double v, double d1,\r
+double calcParIntO2(double b, double d, double t, double v, double d1,\r
double d2, double d3, double eta) {\r
\r
//Elmente anordnen\r
\r
}\r
\r
-double slpR::calcOrtIntO2(double b, double d, double t, double v, double d1,\r
+double calcOrtIntO2(double b, double d, double t, double v, double d1,\r
double d2, double d3, double eta) {\r
\r
//Elmente anordnen\r
if ((b * b + d * d) * eta <= d1 * d1 + d2 * d2 + d3 * d3) {\r
return calcOrtIntQX1X2(b, d, t, v, d1, d2, d3);\r
} else {\r
- cout << "2error";\r
+ // cout << "2error";\r
return calcOrtIntA(b, d, t, v, d1, d2, d3);\r
}\r
}\r
\r
-double slpR::calcParIntO3(double b, double d, double t, double v, double d1,\r
+double calcParIntO3(double b, double d, double t, double v, double d1,\r
double d2, double d3, double eta) {\r
\r
//Achsen vertauschen\r
\r
}\r
\r
-double slpR::calcOrtIntO3(double b, double d, double t, double v, double d1,\r
+double calcOrtIntO3(double b, double d, double t, double v, double d1,\r
double d2, double d3, double eta) {\r
\r
if (max(b, t) * eta > d1) {\r
}\r
}\r
\r
-double slpR::cParO0(double b, double d, double t, double v, double d1,\r
+double cParO0(double b, double d, double t, double v, double d1,\r
double d2, double d3, double* eta) {\r
double tmp = 0;\r
\r
return 0;\r
}\r
\r
-double slpR::cParO1(double b, double d, double t, double v, double d1,\r
+double cParO1(double b, double d, double t, double v, double d1,\r
double d2, double d3, double* eta) {\r
return calcParIntA(b, d, t, v, d1, d2, d3);\r
}\r
\r
-double slpR::cParO2(double b, double d, double t, double v, double d1,\r
+double cParO2(double b, double d, double t, double v, double d1,\r
double d2, double d3, double* eta) {\r
\r
double tmp = 0;\r
return 0;\r
}\r
\r
-double slpR::cParO3(double b, double d, double t, double v, double d1,\r
+double cParO3(double b, double d, double t, double v, double d1,\r
double d2, double d3, double* eta) {\r
\r
double tmp = 0;\r
#ifndef HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
#define HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
-namespace slpR {
+//namespace slpR {
// sol = g0(p,y,x,l);
// double g_AY(double, double, double, double);
double*);
-}
+//}
#endif
% Test ausführen
%Anzahl der Schritte
-steps = 11;
+steps = 40;
%LShape adaptiv anisotrop
-A_run('exmpl_2DLShape', steps, 0, 1, 0.5, 0.5, 'testAA_')
+A_run('exmpl_2DQuad', steps, 0.7, [2 1], 0.5, 0.5, 'testAA_')
%LShape adaptiv isotrop
% A_run('exmpl_2DLShape', steps, 0, 0.5, 0, 0, 'testAI_')
%Ergebnis Plotten und Teil der Daten ausgeben sowie abspeichern der figure
-A_plots({['meshSave/testAA_1_' num2str(steps)]},'figure')
+A_plots({['meshSave/testAA_1_' num2str(steps)]},'exmplAA_2DQuad')
projects / bacc.git / commitdiff