From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Sat, 21 May 2011 16:55:57 +0000 (+0000)
Subject: [src] Matlab g G und F entfernt
X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=ca6c2a19d5658cb76528354c62d40749b3a91fe8;p=bacc.git

[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
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
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
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 49175fd..4bf70d2 100644
--- 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)
 		//		d[i] = y0[i]-x0[i];
 
-			printf("(%d,%d)",rx,ry);
+			//printf("(%d,%d)",rx,ry);
 			if (rx == ry) {
-				printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
-				printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);
-				printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);
-				printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
+				//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
+				//printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);
+				//printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);
+				//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
 				if (rxa == rya) {
 					//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],
 					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
diff --git a/src/exmpl_2DLShape.mat b/src/exmpl_2DLShape.mat
index feb1a94..655e7bd 100644
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 eec4281..994a162 100644
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 a00867a..d2e8234 100644
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
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 bef2205..8502634 100644
--- 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
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 4eca864..9730b85 100644
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -1,51 +1,62 @@
 
+%% Laden des Gitters
+load exmpl_2DLShape %Energienorm sollte gegen 8.28466 konvergieren
 
-load exmpl_2DLShape
-
+%% Rotationen einzelnen von Elementen (zum Testen)
 % coordinates = coordinates(:,[ 1 3 2]);
-elements(2,:) = elements(2,[ 3 4 1 2]);
-elements(3,:) = elements(3,[ 3 2 1 4]);
+% elements(2,:) = elements(2,[ 3 4 1 2]);
+% elements(3,:) = elements(3,[ 3 2 1 4]);
 
+%% Verfeinerungen (isotrop)
+for i = 1:0
+  [coordinates,elements]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
+end
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
-% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
-% build_A(coordinates,elements);
 
-% A = build_A2(coordinates,elements);
-% sum(sum (A-A'));
-% b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));
-% x = A\b;
-% xe = x'*A*x
+for ref = 1:4
 
-A = build_A(coordinates,elements)
-sum(sum (A-A'));
-b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));
-x = A\b;
-xe = x'*A*x
+  %% Testen mit Quadratur... und mex_G00
+  % A = build_A2(coordinates,elements);
+  % sum(sum (A-A'));
+  % b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));
+  % x = A\b;
+  % xe = x'*A*x
 
-% [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
-% A_fine = build_A2(coordinates_fine,elements_fine)
-% sum(sum(A_fine-A_fine'))
-% b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2));
-% x_fine = A_fine\b;
-% xe_fine = x_fine'*A_fine*x_fine;
-% 
-% x_interpol = zeros(size(elements_fine,1),1);
-% for i = 1:size(elements,1)
-%     x_interpol(f2s(i,:)) = x(i);
-% end
-% xe_interpol = x_interpol'*A_fine*x_interpol;
-% 
-% [x_fine x_interpol x_fine-x_interpol]
-% 
-% % l = x_fine(f2s(1,:));
-% % T4 = [1 1 1 1; 1 -1 1 -1; 1 1 -1 -1;1 -1 -1 1];
-% % l
-% % T4*l*1/4
-% 
-% [xe xe_fine xe_interpol]
+  %% Testen mit MEX -> build_A
+  A = build_A(coordinates,elements);
+  sum(sum (A-A'));
+  b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));
+  x = A\b;
+  xe = x'*A*x;
+
+  %% isotrope Verfeinerung und Aufbauen
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
+  A_fine = build_A(coordinates_fine,elements_fine);
+  sum(sum(A_fine-A_fine'));
+  b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2));
+  x_fine = A_fine\b;
+  xe_fine = x_fine'*A_fine*x_fine;
+
+
+  x_interpol = zeros(size(elements_fine,1),1);
+  for i = 1:size(elements,1)
+      x_interpol(f2s(i,:)) = x(i);
+  end
+  xe_interpol = x_interpol'*A_fine*x_interpol;
+
+  % [x_fine x_interpol x_fine-x_interpol]
+
+  % Wie soll verfeinert werden
+  marked = refineType(x_fine(f2s)',0.3);
+
+  [xe xe_fine xe_interpol]
+
+  [coordinates, elements] = refineQuad(coordinates,elements,marked);
+
+end
 
 % re_times = 6;
 % r_type = 1;
@@ -68,4 +79,3 @@ xe = x'*A*x
 % legend('x^TAx')
 % xlabel('Elemente')
 % ylabel('xNorm')
-