+++ /dev/null
-function A_setParam(paras)
-%setzt die Parameter fuer die Durchlaeufe
-
-global G_para;
-G_para = paras;
-
-end
-
+++ /dev/null
-function [coo, ele] = exportCOEL(coo, ele, file)
-%EXPORTCOEL Summary of this function goes here
-% Detailed explanation goes here
-
-if(size(ele,2) == 3)
- disp 'erweitere auf 4 Koordinaten pro Element'
- ele = ele(:,[1 2 3 3]);
-
- for i = 1:size(ele,1);
-
- tri = ele(i,:);
- a = (coo(tri(2),:)-coo(tri(1),:));
- x = coo(tri(3),:) + a;
-
- found = 0;
- for j = 1:size(coo,1)
- if(coo(j,:) == x)
- found = j;
- break;
- end
- end
-
- if(found == 0)
- coo(j+1,:) = x;
- ele(i,3) = j+1;
- else
- ele(i,3) = j;
- end
- end
-end
-if(size(ele,2)==4)
- disp 'schreibe in Datei'
- dat = fopen(file,'w');
-
- for i=1:size(coo,1);
- fprintf(dat,'v %.4f %.4f %.4f\n',coo(i,:));
- end
- fprintf(dat,'\n');
- for i=1:size(ele,1);
- fprintf(dat,'f %.0f %.0f %.0f %.0f\n',ele(i,:));
- end
-
- fclose(dat);
-end
-end
-
+++ /dev/null
-function [nodes,weights] = gauss(n,varargin)
-
-% Compute Gauss-Legendre quadrature on compact interval [a,b]
-%
-% Usage: [NODES,WEIGHTS] = gauss(N [,A,B])
-%
-% Here, N is the length of the Gaussian quadrature rule, and the
-% optional scalars A and B determine the compact interval. By default,
-% there holds A = -1 and b = +1.
-%
-% The function returns an (N x 1) column vector WEIGHTS containing the
-% quadrature weights and an (1 x N) row vector NODES containing the
-% corresponding quadrature nodes.
-%
-% Example: Suppose we aim to approximate the integral
-%
-% int_a^b f dx
-%
-% Assume that F is a Matlab function, which takes a column vector of
-% evaluation points X and returns a column vector Y of function values,
-% where Y(j) = F(X(j)). Then, the numerical integration reads
-%
-% [NODES,WEIGHTS] = gauss(n,A,B);
-% INT = WEIGHTS*F(NODES);
-%
-% Author: Dirk Praetorius - 11.03.2010
-
-
-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;
-
-if nargin >= 3
- a = varargin{1};
- b = varargin{2};
- weights = 0.5*abs(b-a)*weights;
- nodes = 0.5*( a+b + nodes*(b-a) );
-end
+++ /dev/null
-function str = gaussgen(start,stop,steps)
-
-str =['{' char(10)];
-for i = steps(1:end)
- str = [str '{' char(10)];
- [n w] = gauss(i,start,stop);
- for j = 1:length(n)
- str = [str char(9) '{' num2str(n(j),16) ',' num2str(w(j),16) '},' char(10)];
- end
- str = [str(1:end-2) '},' char(10)];
-end
-
-
-
-str = [str(1:end-2) '}'];
-end
\ No newline at end of file
+++ /dev/null
-function neigh = genNeig(ele)
-%momentan nur fuer Netze ohne HangingNodes
-
-neigh = zeros(size(ele,1),8);
-
-
- for i = 1:size(ele,1)
-
- list = [];
- n = 1;
- for k = 1:4
- [r c] = find(ele==ele(i,k));
- list = [list r'];
- end
- list(list~=i)
- list = sort(list(list~=i));
-
- li = [];
- for k = 1:length(list)-1
- if(list(k)==list(k+1))
- neigh(i,n) =list(k);
- n = n+1;
- k = k+1;
- else
- li = [li list(k)];
- end
- end
-
- % li
- % neigh(i,:)
- end
-end
-
\ No newline at end of file
+++ /dev/null
-/*\r
- * Art der Berechnung durch Parameter bestimmt...\r
- * Analytisch/Semianalytisch/...\r
- *\r
- *\r
- */\r
-\r
-#include <iostream>\r
-#include <cmath>\r
-#include <cassert>\r
-#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
-\r
-#include "slpRectangle.hpp"\r
-\r
-using namespace std;\r
-using namespace slpR;\r
-\r
-int dimOfVec(double* vec) {\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
- return -1;\r
- }\r
-\r
-}\r
-\r
-inline int dimOfThird(int a, int b) {\r
- return ((-(a + b) % 3) + 3) % 3;\r
-}\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
-\r
- //Sicherheitsabfragen zu Datengroessen\r
- if ((nrhs != 4))\r
- mexErrMsgTxt(\r
- "expected (coordinates(Nx3),elements(Mx4),mu(double),type(int))");\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
- if (cn != 3)\r
- mexErrMsgTxt("expected coordinates (Nx3)");\r
- cout << " Coordinaten:" << 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
-\r
- //Auslesen der Parameter\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 mu = *(mxGetPr(prhs[2]));\r
- int type = (int) *(mxGetPr(prhs[3]));\r
-\r
- //Initialisieren der Hilfsvektoren\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 * 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 * 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
- // Welche Funktion soll verwendet werden\r
- double(*ctypeP)(double, double, double, double, double, double, double,\r
- double);\r
- double(*ctypeO)(double, double, double, double, double, double, double,\r
- double);\r
-\r
- //Art der Berechnung bestimmen\r
- cout << " Type:" << type << endl;\r
- switch (type) {\r
- case 0:\r
- ctypeP = calcParIntO0;\r
- ctypeO = calcOrtIntO0;\r
- break;\r
- case 1:\r
- ctypeP = calcParIntO1;\r
- ctypeO = calcOrtIntO1;\r
- break;\r
- case 2:\r
- ctypeP = calcParIntO2;\r
- ctypeO = calcOrtIntO2;\r
- break;\r
- case 3:\r
- ctypeP = calcParIntO3;\r
- ctypeO = calcOrtIntO3;\r
- break;\r
-\r
- }\r
-\r
- //LageInformationen\r
- int rx, rxa, rxb, ry, rya, ryb;\r
-\r
-// cout << " Progress: #";\r
- count = 0;\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
- // Lageeigenschaften des Flaechenstuecks\r
-\r
- rxa = dimOfVec(xa);\r
- rxb = dimOfVec(xb);\r
- rx = dimOfThird(rxa, rxb);\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
- }\r
- }\r
-\r
- for (k = i+1; k < em; ++k) {\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
- // Lageeigenschaften des Flaechenstuecks\r
-\r
- rya = dimOfVec(ya);\r
- ryb = dimOfVec(yb);\r
- ry = dimOfThird(rya, ryb);\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
- }\r
- }\r
-\r
- if (rx == ry) {\r
- if (rxa == rya) {\r
- tmp\r
- = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]),\r
- fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],\r
- d[rxb], d[rx], mu);\r
-\r
- } else {\r
- tmp\r
- = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]),\r
- fabs(yb[rxa]), fabs(ya[rxb]), d[rxa],\r
- d[rxb], d[rx], mu);\r
- }\r
-\r
- } else {\r
-\r
- if (rxa == rya) {\r
- tmp\r
- = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]),\r
- fabs(ya[rya]), fabs(yb[ryb]), d[rxb],\r
- d[rxa], d[rx], mu);\r
- } else if (rxa == ryb) {\r
- tmp\r
- = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]),\r
- fabs(yb[ryb]), fabs(ya[rya]), d[rxb],\r
- d[rxa], d[rx], mu);\r
- } else if (rxb == rya) {\r
- tmp\r
- = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]),\r
- fabs(ya[rya]), fabs(yb[ryb]), d[rxa],\r
- d[rxb], d[rx], mu);\r
- } else {\r
- tmp\r
- = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]),\r
- fabs(yb[ryb]), fabs(ya[rya]), d[rxa],\r
- d[rxb], d[rx], mu);\r
- }\r
-\r
- }\r
- A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
- // if(k!=j)\r
-// A[(j * em) + k] = 1. / (4 * M_PI) * tmp;\r
-/* if(count++ > ((em*(em+1))/2)/10){\r
- count = 0;\r
- cout << "#";\r
- cout.flush();\r
- }*/\r
-\r
- }\r
-\r
- }\r
- cout << endl;\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(d);\r
- //mxFree(dt);\r
-\r
- return;\r
-}\r
+++ /dev/null
-#include <iostream>
-#include "mex.h"
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
- std::cout << "Hello" << std::endl;
-
-
-
-}
+++ /dev/null
-function [ y ] = neville2(x,y,t )\r
-%[ p ] = neville2(x,y,t )\r
-\r
-\r
- n=length(y);\r
- t(2:n) = t(1);\r
- \r
- for m = 2:n\r
- \r
- j=n-m+1;\r
- %... = zeigt MATLAB das die Zeile noch weiter geht\r
- y = ((t(1:j) - x(1:j)) .* y(2:j+1) ... \r
- -(t(1:j) - x(m:n)) .* y(1:j)) ...\r
- ./ (x(m:n)-x(1:j));\r
- \r
- end\r
-\r
-\r
-end
\ No newline at end of file
+++ /dev/null
-function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
-%
-% in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
-%
-% Berechnet das Vierfachintegral über die Funktion f(x1,x2,y1,y2) mit den
-% Grenzen [a b]x[c d]x[r t]x[u v]. Dazu wird eine Gaussquadratur über w
-% Punkte verwendet.
-%
-% P. Schaefer
-
- [x1n x1w] = gauss(w,a,b);
- [x2n x2w] = gauss(w,c,d);
- [y1n y1w] = gauss(w,r,t);
- [y2n y2w] = gauss(w,u,v);
-
- in = 0;
- for i=1:length(y2n)
- for j=1:length(y1n)
- for k=1:length(x2n)
- for l = 1:length(x1n)
-
- in = in + y2w(i) * y1w(j) * x2w(k) * x1w(l) *...
- f(x1n(l),x2n(k),y1n(j),y2n(i));
-
- end
- end
- end
- end
-
-end
\ No newline at end of file
+++ /dev/null
-function in = surfQuad(f,a,b,c,d,v)\r
-%\r
-% in = surfQuad(f,a,b,c,d,v)\r
-%\r
-% Berechnet das Doppelintegral über die Funktion f(x,y) mit den Grenzen\r
-% [a b]x[c d]. Dazu wird eine Gaussquadratur über v Punkte verwendet.\r
-%\r
-% P. Schaefer\r
-\r
- [xn xw] = gauss(v,a,b);\r
- [yn yw] = gauss(v,c,d);\r
- \r
- in = 0;\r
- for i = 1:length(yn)\r
- for j = 1:length(xn)\r
- in = in + yw(i) * xw(j)*f(xn(j),yn(i));\r
- end\r
- end\r
-\r
-end
\ No newline at end of file
+++ /dev/null
-\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 Uebersicht\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 A3 = mex_build_AU(coordinates,elements,1,3);\r I\r dat(I,1:4) = [h A0(1,2) A2(1,2) A3(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','Quad Achse','location','best');\rxlabel('Elementgroesse (k??rzeste Seite)');\rylabel('Integral');\rtitle('Integral bei kleiner werdenden Element');\r\r%% Netzverfeinerung mit schlechtem Netz\r\r% h = (1/2)^25;\r% coordinates=coo(h,diff);\r% \r% A_loadMesh(coordinates,elements,neigh);\r% datA=[];\r% for I = 1:4\r% datA(I,:) = A_stepIso(1,[0 3 2]);\r% % datA(I,:) = A_stepAniso(1,[0 3 2],1,0);\r% I\r% end\r% figure(3)\r% loglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,9),datA(2:end,1),datA(2:end,6))\r% \r% legend('Analytisch','Quad Element','Quad Achse','location','best');\r% xlabel('Anzahl der Elemente');\r% ylabel('mu Schaetzer');\r% title('mu Schaetzer mit "schlechtem" Startnetz');\r% \r% figure(4)\r% loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7))\r% \r% legend('Analytisch','Quad Element','Quad Achse','location','best');\r% xlabel('Anzahl der Elemente');\r% ylabel('EnergieNorm ^2 ');\r% title('EnergieNorm ^2 mit "schlechtem" Startnetz');\r% \r% datA\r
\ No newline at end of file
+++ /dev/null
-
-
-coo = @(h,diff)[0 0 0;1 0 0; 1 1 0;0 1 0; 0 0 0 ; 0+h 0 0; 0+h 1 0 ; 0 1 0]+[zeros(4,3);repmat(diff,4,1)];
-elements=[1 2 3 4;5 6 7 8];
-neigh = zeros(2,8);
-dat = [];
-
-diff = [5 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);
- 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;
-end
-
-figure(2)
-loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),'-.',dat(:,1),abs(dat(:,4)),'--');
-
-legend('Analytisch','Quad Element','Quad Achse','location','best');
-xlabel('Elementgroesse (k??rzeste Seite)');
-ylabel('Integral');
-title('Integral bei kleiner werdenden Element');
-
-%% Netzverfeinerung mit schlechtem Netz
-
-% h = (1/2)^25;
-% coordinates=coo(h,diff);
-%
-% A_loadMesh(coordinates,elements,neigh);
-% datA=[];
-% for I = 1:4
-% datA(I,:) = A_stepIso(1,[0 3 2]);
-% % datA(I,:) = A_stepAniso(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,9),datA(2:end,1),datA(2:end,6))
-%
-% legend('Analytisch','Quad Element','Quad Achse','location','best');
-% xlabel('Anzahl der Elemente');
-% ylabel('mu Schaetzer');
-% title('mu Schaetzer mit "schlechtem" Startnetz');
-%
-% figure(4)
-% loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7))
-%
-% legend('Analytisch','Quad Element','Quad Achse','location','best');
-% xlabel('Anzahl der Elemente');
-% ylabel('EnergieNorm ^2 ');
-% title('EnergieNorm ^2 mit "schlechtem" Startnetz');
-%
-% datA
+++ /dev/null
-
-
-coo = @(h,diff)[0 0 0;h 0 0; h 1 0;0 1 0; 0 0 0 ; 0+h 0 0; 0+h 1 0 ; 0 1 0]+[zeros(4,3);repmat(diff,4,1)];
-elements=[1 2 3 4;5 6 7 8];
-neigh = zeros(2,8);
-dat = [];
-
-diff = [5 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(1,[1:4,1])',:);
-fill3(current(:,1),current(:,2),current(:,3),'r');
-current = coordinates(elements(2,[1:4,1])',:);
-fill3(current(:,1),current(:,2),current(:,3),'y');
-hold off
-legend('Element1(h)','Element2(h)','Element1(h/2)','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);
- 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;
-end
-
-figure(2)
-loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),'-.',dat(:,1),abs(dat(:,4)),'--');
-
-legend('Analytisch','Quad Element','Quad Achse','location','best');
-xlabel('Elementgroesse (k??rzeste Seite)');
-ylabel('Integral');
-title('Integral bei kleiner werdenden Element');
-
-%% Netzverfeinerung mit schlechtem Netz
-%
-% h = (1/2)^25;
-% coordinates=coo(h,diff);
-%
-% A_loadMesh(coordinates,elements,neigh);
-% datA=[];
-% for I = 1:4
-% datA(I,:) = A_stepIso(1,[0 3 2]);
-% % datA(I,:) = A_stepAniso(1,[0 3 2],1,0);
-% I
-% end
-% figure(3)
-% loglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,9),datA(2:end,1),datA(2:end,6))
-%
-% legend('Analytisch','Quad Element','Quad Achse','location','best');
-% xlabel('Anzahl der Elemente');
-% ylabel('mu Schaetzer');
-% title('mu Schaetzer mit "schlechtem" Startnetz');
-%
-% figure(4)
-% loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7))
-%
-% legend('Analytisch','Quad Element','Quad Achse','location','best');
-% xlabel('Anzahl der Elemente');
-% ylabel('EnergieNorm ^2 ');
-% title('EnergieNorm ^2 mit "schlechtem" Startnetz');
-%
-% datA
+++ /dev/null
-
-
-elements=[1 2 3 4;5 6 7 8];
-neigh = zeros(2,8);
-
-dat = [];
-
-h = 1/2;
-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];
-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 0 h; 2 h h ; 2 h 0];
-current = coordinates(elements(2,[1:4,1])',:);
-fill3(current(:,1),current(:,2),current(:,3),'r');
-hold off
-legend('Element1(h/2)','Element2(h/2)','Element1(h/4)','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 0 h; 2 h h ; 2 h 0];
-
- A0 = mex_build_AU(coordinates,elements,0,0);
- A1 = mex_build_AU(coordinates,elements,1,1);
- A2 = mex_build_AU(coordinates,elements,1,2);
- I
- dat(I,1:4) = [h A0(1,2) A1(1,2) A2(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','Element vertauschen','Quad Element','location','southeast');
-
-xlabel('Elementgroesse');
-ylabel('Integral');
\ No newline at end of file
+++ /dev/null
-figure(1)
-plotShape(coordinates,elements);
-plotMark(1,coordinates,elements);
-
-[coordinates elements neigh]= refineQuad(coordinates,elements,neigh,[3 ones(1,size(elements,1)-1)]);
-
-figure(2)
-plotShape(coordinates,elements);
+++ /dev/null
-
-N = size(elements,1);
-
-tic
-A = mex_build_A(coordinates,elements);
-t = toc;
-b = sqrt(sum(quadNorm(coordinates,elements,'w').^2,2));
-
-x = A\b;
-xe = x'*A*x;
-
-x;
-xe
-t
\ No newline at end of file
+++ /dev/null
-\r
-dataIso =[];\r
-dataAniso =[];\r
-maxSize = 10^2;\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
-%% Isotrope Verfeinerung\r
-disp 'Isotrop'\r
-load(Netz)\r
-\r
-ref = 0;\r
-while size(elements,1)<maxSize\r
- tic\r
- ref = ref+1;\r
-\r
- [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
- A_fine = mex_build_A(coordinates_fine,elements_fine);\r
- b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
- x_fine = A_fine\b;\r
- xe_fine = x_fine'*A_fine*x_fine;\r
-\r
- ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
- marked = 2*ones(1,size(elements,1));\r
-\r
- dataIso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
- [coordinates, elements] = refineQuad(coordinates,elements,marked);\r
- toc\r
-end\r
-dataIso\r
-\r
-\r
-%% Anisotrope Verfeinerung\r
-disp 'Anisotrop'\r
-load(Netz)\r
-\r
-ref = 0;\r
-while size(elements,1)<maxSize\r
- tic\r
- ref = ref+1;\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
- toc\r
-end\r
-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(dataIso(:,1),dataIso(:,2),dataAniso(:,1),dataAniso(:,2));\r
-\r
-legend('Isotrop','Anisotrop')\r
-xlabel('Elemente')\r
-ylabel('Fehler')\r
-\r
-figure (2)\r
-loglog(dataIso(:,1),dataIso(:,3),dataAniso(:,1),dataAniso(:,3));\r
-\r
-legend('Isotrop','Anisotrop')\r
-xlabel('Elemente')\r
-ylabel('EnergieNorm')\r
+++ /dev/null
-\r
-dataIso =[];\r
-dataAniso =[];\r
-maxSize = 500;\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
-%% Anisotrope Verfeinerung\r
-disp 'Anisotrop'\r
-load(Netz)\r
-\r
-ref = 0;\r
-while size(elements,1)<maxSize\r
- tic\r
- ref = ref+1; \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
- save ([file], 'coordinates', 'elements')\r
- toc\r
-end\r
-dataAniso\r
-\r
-clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine\r
-\r
+++ /dev/null
-\r
-mex mex_build_As1.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
-if(~exist('d1'))\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_As1(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
+++ /dev/null
-\r
-mex mex_build_AU.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 16.2265 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
- %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_AU(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
+++ /dev/null
-function r = testm(a,b)
-
-global A;
-A = a;
-te(1)
-
-r = A;
-
-end
-
-
-function b = te(b)
-
-global A;
-
-disp(A);
-
-A = A+1;
-
-end
\ No newline at end of file
+++ /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