[src] Matlab g G und F entfernt

[src] Art der Verfeinerung Implementiert
[src] refine an 4 Koo pro Ele angepasst
[src] test_solveE adaptive Verfeinerung (ob verfeinert wird oder nicht... muss noch entschieden werden)
[src] Beispielgitter in 4Koo pro Ele übersetzt mit expCOEL

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

diff --git a/src/F_ort.m b/src/F_ort.m
deleted file mode 100644 (file)
index 2c5636b..0000000
--- a/src/F_ort.m
+++ /dev/null
@@ -1,14 +0,0 @@
-function sol = F_ort(x, y, delta)
-
-    xd = x(1) -delta(1);
-    xyd=x(2)-y(1)-delta(2);
-    yd = y(2)+delta(3);
-    
-    sol = -G00(0.5,[y(2) x(1)]',[-delta(3) delta(1)],xyd)...
-        -xd*xyd*G00(-0.5,[x(2) y(2)],[yd -delta(3)],xd)...
-        +xd*g0(0.5,y(2),-delta(3),norm([xd xyd]))...
-        -yd*xyd*G00(-0.5,[x(1) x(2)],[delta(1) yd],-yd)...
-        +yd*g0(0.5,x(1),delta(1),norm([xyd yd]));
-
-end
-

diff --git a/src/F_par.m b/src/F_par.m
deleted file mode 100644 (file)
index c84bb42..0000000
--- a/src/F_par.m
+++ /dev/null
@@ -1,22 +0,0 @@
-function sol = F_par(x1,x2,y1,y2,d1,d2,d3)
-%     fprintf('%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f',x1,x2,y1,y2,d1,d2,d3)
-%     xyd=x-y-delta(1:2);
-    
-    sol = (x1-y1-d1)*(x2-y2-d2);
-    
-    if(sol~=0)
-        sol = sol * G00(-0.5,x1,x2,y1+d1,y2+d2,d3);
-    end
-    
-    if((x1-y1-d1)~=0)
-        sol = sol - (x1-y1-d1)*g0(0.5,x1,y1+d1,sqrt((x2-y2-d2)^2+d3^2));
-    end
-    
-    if((x2-y2-d2)~=0)
-        sol = sol -(x2-y2-d2)*g0(0.5,x2,y2+d2,sqrt((x1-y1-d1)^2+d3^2));
-    end
-    
-    sol = sol + 1/3*((x1-y1-d1)^2+(x2-y2-d2)^2+d3^2)^(1.5);
-
-end
-

diff --git a/src/G00.m b/src/G00.m
deleted file mode 100644 (file)
index b725d43..0000000
--- a/src/G00.m
+++ /dev/null
@@ -1,26 +0,0 @@
-function sol = G00(p,y1,y2,x1,x2,l)
-% fprintf('%.1f | %.1f %.1f | %.1f %.1f | %.1f +',p,x1,x2,y1,y2,l);
-   sol = 0;
-    if(p==-1.5)
-       if(l==0)
-          sol = -sqrt((y1-x1)^2+(y2-x2)^2)/((y1-x1)*(y2-x2));
-       else
-          sol = sign((y1-x1)*(y2-x2))/(2*abs(l)) ...
-              * acos(-2*(y1-x1)^2*(y2-x2)^2 ...
-              /(((y1-x1)^2+l^2)*((y2-x2)^2+l^2))+1);
-       end
-    elseif(p==-0.5)
-       if(l~=0)
-          sol = 2*p*l^2*G00(p-1,y1,y2,x1,x2,l);
-       end
-       if((y1-x1)~=0)
-          sol = sol + (y1-x1)*g0(p,y2,x2,sqrt((y1-x1)^2+l^2));
-       end
-       if((y2-x2)~=0)
-          sol = sol + (y2-x2)*g0(p,y1,x1,sqrt((y2-x2)^2+l^2));
-       end
-       sol = sol / (2*p+2);
-    else
-        error 'no case for p defined'
-    end
-end

diff --git a/src/build_A.cpp b/src/build_A.cpp
index 49175fdde8e943abce4f8b4c460e0520660a069b..4bf70d2a11b09e20824effbf0c8f4023bcb6291e 100644 (file)
--- a/src/build_A.cpp
+++ b/src/build_A.cpp
@@ -216,12 +216,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                //      for (i = 0; i<3;++i)\r
                //              d[i] = y0[i]-x0[i];\r
 \r
-                       printf("(%d,%d)",rx,ry);\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
+                               //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

diff --git a/src/exmpl_2DLShape.mat b/src/exmpl_2DLShape.mat
index feb1a9435624f9e922a173ddea746e683980b0e6..655e7bde2958eaedb65001b88e38609ef0253cdb 100644 (file)
Binary files a/src/exmpl_2DLShape.mat and b/src/exmpl_2DLShape.mat differ

diff --git a/src/exmpl_2DQuad.mat b/src/exmpl_2DQuad.mat
index eec428172f09e2c930f76b2658f173e25136910e..994a16200fdc7e0fa01c5cd311d3fb83afe68323 100644 (file)
Binary files a/src/exmpl_2DQuad.mat and b/src/exmpl_2DQuad.mat differ

diff --git a/src/exmpl_3DCube.mat b/src/exmpl_3DCube.mat
index a00867af717ca3ed6d408f11ebe9959159d34bd6..d2e82344d2a377bd00f5c8a196cb38fe44ee98e0 100644 (file)
Binary files a/src/exmpl_3DCube.mat and b/src/exmpl_3DCube.mat differ

diff --git a/src/g0.m b/src/g0.m
deleted file mode 100644 (file)
index 023fa2d..0000000
--- a/src/g0.m
+++ /dev/null
@@ -1,34 +0,0 @@
-function sol = g0(p, y, x, l)
- %fprintf('%.1f | %.1f | %.1f | %.1f +',p,x,y,l);
- %sol = nan;
-     yx = y-x;
-     if(l~=0)   
-        if(p==0.5)
-            sol = yx/2*((yx^2 + l^2)^0.5)+...
-                l^2/2*arsinh(yx/abs(l));
-      %      fprintf('%.2f |',sol);
-        elseif(p==0)
-            sol = yx;
-        elseif(p==-0.5)
-            sol = arsinh(yx/abs(l));
-        elseif(p==-1)
-            sol = atan(yx/abs(l));
-        elseif(p==-1.5)
-            sol =  yx/(l^2)*((yx^2 + l^2)^-0.5);
-        else
-            sol = yx*((yx^2 + l^2)^p)+2*p*l^2*g0(p-1,y,x,l)/(2*p+1);
-        end
-     else
-        if(p==-0.5)
-            sol = sign(yx)*log(abs(yx));
-        else
-            sol = (yx).*abs(yx).^(2*p)/(2*p+1);
-        end
-     end
-    
-end
-
-function sol = arsinh(x)    %inline
-    sol = log(x+norm([x 1]));
-end
-

diff --git a/src/refineQuad.m b/src/refineQuad.m
index bef2205bcfa3c9d4f36338d24b36c57cb7b7002d..85026348552f7cfb756a0bd0a7b5579c03015dcc 100644 (file)
--- a/src/refineQuad.m
+++ b/src/refineQuad.m
@@ -14,47 +14,39 @@ for i = 1:c_loop
     
     
     el = elements(i,:);
-    if(type(i)==0)
+    if(type(i)==1)
         continue;
-    elseif(type(i)==1)
-        newc1 = (coordinates(el(2),:)+coordinates(el(1),:))/2;
-        newc2 = (coordinates(el(3),:)+coordinates(el(1),:))/2;
-        coordinates(c_coo+1,:) = newc1;
-        coordinates(c_coo+2,:) = newc2;
-        elements(i,:) = [el(1),c_coo+1,c_coo+2];
-        
-        newc3 = (coordinates(el(2),:)+coordinates(el(3),:))/2;
-        coordinates(c_coo+3,:) = newc3;
-        elements(c_ele+1,:) = [c_coo+1,el(2),c_coo+3];
-        elements(c_ele+2,:) = [c_coo+2,c_coo+3,el(3)];
-        
-        newc4 = newc2+(coordinates(el(2),:)-coordinates(el(1),:));
-        newc5 = newc1+(coordinates(el(3),:)-coordinates(el(1),:));
-        coordinates(c_coo+4,:) = newc4;
-        coordinates(c_coo+5,:) = newc5;
-        elements(c_ele+3,:) = [c_coo+3,c_coo+4,c_coo+5];
+    elseif(type(i)==2)
+        coordinates(c_coo+1,:) = (coordinates(el(1),:)+coordinates(el(2),:))/2;
+        coordinates(c_coo+2,:) = (coordinates(el(2),:)+coordinates(el(3),:))/2;
+        coordinates(c_coo+3,:) = (coordinates(el(3),:)+coordinates(el(4),:))/2;
+        coordinates(c_coo+4,:) = (coordinates(el(1),:)+coordinates(el(4),:))/2;
+        coordinates(c_coo+5,:) = (coordinates(el(1),:)+coordinates(el(3),:))/2;
+
+        elements(i,:) = [c_coo+4,c_coo+5,c_coo+3,el(4)];
+        elements(c_ele+1,:) = [c_coo+5,c_coo+2,el(3),c_coo+3];
+        elements(c_ele+2,:) = [c_coo+1,el(2),c_coo+2,c_coo+5];
+        elements(c_ele+3,:) = [el(1),c_coo+1,c_coo+5,c_coo+4];
     
         fa2so(i,2:4)=c_ele+1:c_ele+3;
-    elseif(type(i)==2)
-        newc = (coordinates(el(3),:)+coordinates(el(1),:))/2;
-        coordinates(c_coo+1,:) = newc;
-        elements(i,3) = c_coo+1;
+    elseif(type(i)==3)
+        coordinates(c_coo+1,:) = (coordinates(el(1),:)+coordinates(el(4),:))/2;
+        coordinates(c_coo+2,:) = (coordinates(el(2),:)+coordinates(el(3),:))/2;
         
-        newc = newc+(coordinates(el(2),:)-coordinates(el(1),:));
-        coordinates(c_coo+2,:) = newc;
+        elements(i,1) = c_coo+1;
+        elements(i,2) = c_coo+2;
         
-        elements(c_ele+1,:) = [c_coo+1,c_coo+2,el(3)];       
-        fa2so(i,2)=c_ele+1;
-    elseif(type(i)==3)
-        newc = (coordinates(el(3),:)+coordinates(el(1),:))/2;
-        coordinates(c_coo+1,:) = newc;
-        elements(i,3) = c_coo+1;
+        elements(c_ele+1,:) = [el(1),el(2),c_coo+2,c_coo+1];
+        fa2so(i,[3 4])=c_ele+1;
+    elseif(type(i)==4)
+        coordinates(c_coo+1,:) = (coordinates(el(1),:)+coordinates(el(2),:))/2;
+        coordinates(c_coo+2,:) = (coordinates(el(4),:)+coordinates(el(3),:))/2;
         
-        newc = newc+(coordinates(el(2),:)-coordinates(el(1),:));
-        coordinates(c_coo+2,:) = newc;
+        elements(i,2) = c_coo+1;
+        elements(i,3) = c_coo+2;
         
-        elements(c_ele+1,:) = [c_coo+1,c_coo+2,el(3)];
-        fa2so(i,2)=c_ele+1;
+        elements(c_ele+1,:) = [c_coo+1,el(2),el(3),c_coo+2];
+        fa2so(i,[2 4])=c_ele+1;
     end
 end
 end

diff --git a/src/refineType.m b/src/refineType.m
new file mode 100644 (file)
index 0000000..f87832b
--- /dev/null
+++ b/src/refineType.m
@@ -0,0 +1,25 @@
+function REF = refineType(xF2S,tau)
+%REFINETYPE Summary of this function goes here
+%   Detailed explanation goes here
+
+T4 = [1 1 1 1; 1 -1 1 -1; 1 1 -1 -1;1 -1 -1 1]/4;
+
+REF=ones(1,size(xF2S,2));
+t1=zeros(1,size(xF2S,2));
+t3=t1;
+t4=t1;
+
+%% Muss ueberhaupt verfeinert werden
+
+% t1(7) = 1;
+
+
+%% Wie muss verfeinert werden
+% Ct = T4*xF2S;
+% t3 = (tau*abs(Ct(3,:)) >= sqrt(Ct(2,:).^2+Ct(4,:).^2));
+% REF(t3) = 3;
+% t4 = (tau*abs(Ct(4,:)) >= sqrt(Ct(2,:).^2+Ct(3,:).^2));
+% REF(t4) = 4;
+REF(~(t4+t3+t1)) = 2;
+end
+

diff --git a/src/test_solveError.m b/src/test_solveError.m
index 4eca86498f829ce9dd05f2b28f2e8a9161e96fad..9730b8573f613f9ba20c212935c7953951008df0 100644 (file)
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -1,51 +1,62 @@
 \r
+%% Laden des Gitters\r
+load exmpl_2DLShape %Energienorm sollte gegen 8.28466 konvergieren\r
 \r
-load exmpl_2DLShape\r
-\r
+%% Rotationen einzelnen von Elementen (zum Testen)\r
 % coordinates = coordinates(:,[ 1 3 2]);\r
-elements(2,:) = elements(2,[ 3 4 1 2]);\r
-elements(3,:) = elements(3,[ 3 2 1 4]);\r
+% elements(2,:) = elements(2,[ 3 4 1 2]);\r
+% elements(3,:) = elements(3,[ 3 2 1 4]);\r
 \r
+%% Verfeinerungen (isotrop)\r
+for i = 1:0\r
+  [coordinates,elements]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
+end\r
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
-% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
-% build_A(coordinates,elements);\r
 \r
-% A = build_A2(coordinates,elements);\r
-% sum(sum (A-A'));\r
-% b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));\r
-% x = A\b;\r
-% xe = x'*A*x\r
+for ref = 1:4\r
 \r
-A = build_A(coordinates,elements)\r
-sum(sum (A-A'));\r
-b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));\r
-x = A\b;\r
-xe = x'*A*x\r
+  %% Testen mit Quadratur... und mex_G00\r
+  % A = build_A2(coordinates,elements);\r
+  % sum(sum (A-A'));\r
+  % b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));\r
+  % x = A\b;\r
+  % xe = x'*A*x\r
 \r
-% [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
-% A_fine = build_A2(coordinates_fine,elements_fine)\r
-% sum(sum(A_fine-A_fine'))\r
-% b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2));\r
-% x_fine = A_fine\b;\r
-% xe_fine = x_fine'*A_fine*x_fine;\r
-% \r
-% x_interpol = zeros(size(elements_fine,1),1);\r
-% for i = 1:size(elements,1)\r
-%     x_interpol(f2s(i,:)) = x(i);\r
-% end\r
-% xe_interpol = x_interpol'*A_fine*x_interpol;\r
-% \r
-% [x_fine x_interpol x_fine-x_interpol]\r
-% \r
-% % l = x_fine(f2s(1,:));\r
-% % T4 = [1 1 1 1; 1 -1 1 -1; 1 1 -1 -1;1 -1 -1 1];\r
-% % l\r
-% % T4*l*1/4\r
-% \r
-% [xe xe_fine xe_interpol]\r
+  %% Testen mit MEX -> build_A\r
+  A = build_A(coordinates,elements);\r
+  sum(sum (A-A'));\r
+  b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));\r
+  x = A\b;\r
+  xe = x'*A*x;\r
+\r
+  %% isotrope Verfeinerung und Aufbauen\r
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
+  A_fine = build_A(coordinates_fine,elements_fine);\r
+  sum(sum(A_fine-A_fine'));\r
+  b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2));\r
+  x_fine = A_fine\b;\r
+  xe_fine = x_fine'*A_fine*x_fine;\r
+\r
+\r
+  x_interpol = zeros(size(elements_fine,1),1);\r
+  for i = 1:size(elements,1)\r
+      x_interpol(f2s(i,:)) = x(i);\r
+  end\r
+  xe_interpol = x_interpol'*A_fine*x_interpol;\r
+\r
+  % [x_fine x_interpol x_fine-x_interpol]\r
+\r
+  % Wie soll verfeinert werden\r
+  marked = refineType(x_fine(f2s)',0.3);\r
+\r
+  [xe xe_fine xe_interpol]\r
+\r
+  [coordinates, elements] = refineQuad(coordinates,elements,marked);\r
+\r
+end\r
 \r
 % re_times = 6;\r
 % r_type = 1;\r
@@ -68,4 +79,3 @@ xe = x'*A*x
 % legend('x^TAx')\r
 % xlabel('Elemente')\r
 % ylabel('xNorm')\r
-\r