[src] Elementdrehen angefangen -> leider Fehlerhaft

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

diff --git a/src/mex_build_As1.cpp b/src/mex_build_As1.cpp
new file mode 100644 (file)
index 0000000..ec82bee
--- /dev/null
+++ b/src/mex_build_As1.cpp
@@ -0,0 +1,329 @@
+#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
+#define EPS 0.02\r
+\r
+using namespace std;\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;\r
+\r
+       int itmp;\r
+\r
+       double * vtmp;\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
+               if (xn[2] != 0)\r
+                       rx = 2;\r
+               else if (xn[1] != 0)\r
+                       rx = 1;\r
+               else\r
+                       rx = 0;\r
+\r
+               if (xa[2] != 0)\r
+                       rxa = 2;\r
+               else if (xa[1] != 0)\r
+                       rxa = 1;\r
+               else\r
+                       rxa = 0;\r
+\r
+               if (xb[2] != 0)\r
+                       rxb = 2;\r
+               else if (xb[1] != 0)\r
+                       rxb = 1;\r
+               else\r
+                       rxb = 0;\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
+               //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
+\r
+               for (k = 0; 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
+                       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
+                       if (yn[2] != 0)\r
+                               ry = 2;\r
+                       else if (yn[1] != 0)\r
+                               ry = 1;\r
+                       else\r
+                               ry = 0;\r
+\r
+                       if (ya[2] != 0)\r
+                               rya = 2;\r
+                       else if (ya[1] != 0)\r
+                               rya = 1;\r
+                       else\r
+                               rya = 0;\r
+\r
+                       if (yb[2] != 0)\r
+                               ryb = 2;\r
+                       else if (yb[1] != 0)\r
+                               ryb = 1;\r
+                       else\r
+                               ryb = 0;\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
+                       //      for (i = 0; i<3;++i)\r
+                       //              d[i] = y0[i]-x0[i];\r
+\r
+\r
+                       if(ry = rx){\r
+\r
+                               if(xa[rxa]*xa[rxa]+xb[rxb]*xb[rxb]>ya[rya]*ya[rya]+yb[ryb]*yb[ryb]){\r
+\r
+       /*                              printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
+                                       printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);\r
+                                       printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);\r
+                                       printf(" | ");\r
+*/\r
+\r
+                                       printf("(%3.d %3.d)",j,k);\r
+\r
+                                       vtmp = xa; xa = ya; ya = vtmp;\r
+                                       vtmp = xb; xb = yb; yb = vtmp;\r
+\r
+                                       itmp = rxa; rxa = rya; rya = itmp;\r
+                                       itmp = rxb; rxb = ryb; ryb = itmp;\r
+\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]",xa[rxa],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
+                       }\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]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);\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
+                                       tmp\r
+                                                       = apply0Int(F_par, 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
+                               } 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
+                                                       = apply0Int(F_par, 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
+\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
+                                                       = apply0Int(F_ort, 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
+                                                       = apply0Int(F_ort, 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
+                                                       = apply0Int(F_ort, 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
+                                                       = apply0Int(F_ort, 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
+               }\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

diff --git a/src/test_solveError.m b/src/test_solveError.m
index 0316541b8239c87108c54fd3fc525a2632731bad..fa537bfa7385c18568f83a1b56573d4da5afb34f 100644 (file)
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -1,7 +1,7 @@
 \r
 dataIso =[];\r
 dataAniso =[];\r
-maxSize = 1050;\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
@@ -40,22 +40,37 @@ while size(elements,1)<maxSize
   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

diff --git a/src/test_solveErrorS1.m b/src/test_solveErrorS1.m
new file mode 100755 (executable)
index 0000000..4143f9c
--- /dev/null
+++ b/src/test_solveErrorS1.m
@@ -0,0 +1,103 @@
+\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_As1(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_As1(coordinates_fine,elements_fine);\r
+  eh = A_fine-mex_build_A(coordinates_fine,elements_fine);\r
+  figure (3)\r
+  spy(eh)\r
+ \r
+  sum(sum(eh))\r
+  if(eh.^2)\r
+    disp 'hä'\r
+  end\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