[src] g G und beide F in Matlab umgesetzt

[src] kleine korrektur an der refineQuad (TypReihenfolge)

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

diff --git a/src/F_ort.m b/src/F_ort.m
new file mode 100644 (file)
index 0000000..2c5636b
--- /dev/null
+++ b/src/F_ort.m
@@ -0,0 +1,14 @@
+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
new file mode 100644 (file)
index 0000000..c214e4d
--- /dev/null
+++ b/src/F_par.m
@@ -0,0 +1,11 @@
+function sol = F_par(x, y, delta)
+
+    xyd=x-y-delta(1:2);
+
+    sol = xyd(1)*xyd(2)*G(-0.5,x,y+delta(1:2),delta(3)...
+        -xyd(1)*g(0.5,x(1),y(1)+delta(1),norm([xyd(2) delta(3)]))...
+        -xyd(2)*g(0.5,x(2),y(2)+delta(2),norm([xyd(1) delta(3)]))...
+        +1/3*(xyd(1)^2+xyd(2)^2+delta(3)^2)^(1.5));
+
+end
+

diff --git a/src/G00.m b/src/G00.m
new file mode 100644 (file)
index 0000000..777e6f5
--- /dev/null
+++ b/src/G00.m
@@ -0,0 +1,18 @@
+function sol = G00(p,y,x,lambda)
+
+    sol = nan;
+    yx = y-x;
+    if(p==-1.5)
+       if(lambda==0)
+           sol = -norm(yx)/prod(yx);
+       else
+          sol = sign(prod(yx))/(2*abs(lambda))*...
+              acos(-2*prod(yx)^2/prod((yx).^2+lambda^2)+1);
+       end
+    elseif(p==-0.5)
+       sol = (2*p*lambda^2*G(p-1,y,x,lambda)...
+           + sum( yx(1).*g(p,y(2),x(2),norm([xy(1) lambda])))...
+           + sum( yx(2).*g(p,y(1),x(1),norm([xy(2) lambda]))))/(2*p+2);
+    end
+
+end

diff --git a/src/g0.m b/src/g0.m
new file mode 100644 (file)
index 0000000..c8e1632
--- /dev/null
+++ b/src/g0.m
@@ -0,0 +1,31 @@
+function sol = g0(p, y, x, lambda)
+
+ sol = nan;
+     yx = y-x;
+     if(lambda~=0)   
+        if(p==0.5)
+            sol = yx/2*((yx^2 + lambda^2)^0.5)+lambda^2/2*arsinh(yx/abs(lambda));
+        elseif(p==0)
+            sol = yx;
+        elseif(p==-0.5)
+            sol = arsinh(yx/abs(lambda));
+        elseif(p==-1)
+            sol = atan(yx/abs(lambda));
+        elseif(p==-1.5)
+            sol =  yx/(lambda^2)*((yx^2 + lambda^2)^-0.5);
+        else
+            sol = yx*((yx^2 + lambda^2)^p)+2*p*lambda^2*g(p-1,y,x,lambda);
+        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/plotShape.m b/src/plotShape.m
index d7c9a9aa643b72ca0cd2d4ceac4f315180ec1eff..8aedde403c1650f85c90d4bd949706c94d2bd712 100644 (file)
--- a/src/plotShape.m
+++ b/src/plotShape.m
@@ -3,11 +3,12 @@ function plotShape(coordinates, elements, varargin)
 % plotShape(coordinates,elements)
 % plotShape(coordinates,elements,'FLAG')
 %
-% Diese Funktion Zeichnet alle Vierecke mit ausgef llten F l c h e n , wenn als
-% dritter Parameter noch 'f' angegeben wird. Werden nur die Kanten
-% gezeichnet. Wenn 'c' angegeben wird werden auch die Ecken als kleine rote
-% Kreise dargestellt
-% Flag: 'c', 'b' oder 'cb' bzw. 'bc'
+% Diese Funktion Zeichnet alle Vierecke mit ausgefuellten Flaechen.
+% FLAG:
+% c -> Koordinaten als rote Kreise darstellen
+% b -> nur Kanten der Elemente einzeichnen
+% n -> Normen auf ein Element einzeichnen (mit Laenge 1)
+% a -> Normen auf ein Element einzeichnen (Laenge durch Flaecheninhalt)
 %
 % P.Schaefer
 
@@ -25,7 +26,9 @@ elseif(optargin==1)
     if(ismember('c',varargin{1}))
         c = 1;
     end
-    if(ismember('n',varargin{1}))
+    if(ismember('a',varargin{1}))
+        n = 2;
+    elseif(ismember('n',varargin{1}))
         n = 1;
     end
 end
@@ -55,7 +58,11 @@ else
 end
 
 if(n)
-    anorm = quad_norm(coordinates,elements);
+    if(n==2)
+        anorm = quad_norm(coordinates,elements,'w');
+    else
+        anorm = quad_norm(coordinates,elements);
+    end
     for idx = 1:eles
         current = sum(coordinates(elements(idx,[2,3])',:),1)/2;
         current = [current ; current+anorm(idx,:)];

diff --git a/src/quad_norm.m b/src/quad_norm.m
index 8fd0c7805d62b8a15e7591094cdba6cbdd24fba3..6dbfe796e2ec2fa7bda47f39f30ccbaff1f3eb9a 100644 (file)
--- a/src/quad_norm.m
+++ b/src/quad_norm.m
@@ -3,8 +3,9 @@ function n = quad_norm(coordinates, elements,varargin)
 % norm = tri_norm(coordinates, elements)
 % norm = tri_norm(coordinates, elements, 'FLAG')
 %
-% This function builds the norm for every triangle
-%
+% Diese Funktion Berechnet die Orthogonalen mit Laenge 1 über alle Flächen
+% FLAG:
+% w -> Laenge entspricht Flaecheninhalt
 % P.Schaefer
 
 %% Parameterueberpruefung

diff --git a/src/refineQuad.m b/src/refineQuad.m
index 3c199fcb04578a7baed263eaac2494af663b0b6d..fffea50f3ab530a51d4469d494524320d0ceceb7 100644 (file)
--- a/src/refineQuad.m
+++ b/src/refineQuad.m
@@ -11,24 +11,6 @@ for i = 1:c_loop
     if(type(i)==0)
         continue;
     elseif(type(i)==1)
-        newc = (coordinates(el(3),:)+coordinates(el(1),:))/2;
-        coordinates(c_coo+1,:) = newc;
-        elements(i,3) = c_coo+1;
-        
-        newc = newc+(coordinates(el(2),:)-coordinates(el(1),:));
-        coordinates(c_coo+2,:) = newc;
-        
-        elements(c_ele+1,:) = [c_coo+1,c_coo+2,el(3)];
-    elseif(type(i)==2)
-        newc = (coordinates(el(3),:)+coordinates(el(1),:))/2;
-        coordinates(c_coo+1,:) = newc;
-        elements(i,3) = c_coo+1;
-        
-        newc = newc+(coordinates(el(2),:)-coordinates(el(1),:));
-        coordinates(c_coo+2,:) = newc;
-        
-        elements(c_ele+1,:) = [c_coo+1,c_coo+2,el(3)];       
-    elseif(type(i)==3)
         newc1 = (coordinates(el(2),:)+coordinates(el(1),:))/2;
         newc2 = (coordinates(el(3),:)+coordinates(el(1),:))/2;
         coordinates(c_coo+1,:) = newc1;
@@ -45,6 +27,24 @@ for i = 1:c_loop
         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)
+        newc = (coordinates(el(3),:)+coordinates(el(1),:))/2;
+        coordinates(c_coo+1,:) = newc;
+        elements(i,3) = c_coo+1;
+        
+        newc = newc+(coordinates(el(2),:)-coordinates(el(1),:));
+        coordinates(c_coo+2,:) = newc;
+        
+        elements(c_ele+1,:) = [c_coo+1,c_coo+2,el(3)];       
+    elseif(type(i)==3)
+        newc = (coordinates(el(3),:)+coordinates(el(1),:))/2;
+        coordinates(c_coo+1,:) = newc;
+        elements(i,3) = c_coo+1;
+        
+        newc = newc+(coordinates(el(2),:)-coordinates(el(1),:));
+        coordinates(c_coo+2,:) = newc;
+        
+        elements(c_ele+1,:) = [c_coo+1,c_coo+2,el(3)];
     end
 end
 end

diff --git a/src/test_refine.m b/src/test_refine.m
index 0446fa06e6552d25d141a2bb2f75e379d9b947b5..5fd7113aa548970c954f0a63d09c08aa01f4054a 100644 (file)
--- a/src/test_refine.m
+++ b/src/test_refine.m
@@ -11,7 +11,7 @@
  
     if(r_type > 3)
         marked = round(3*rand(1,eles));
-    elseif(r_type < 0)
+    elseif(r_type <= 0)
         return;
     else
         marked = ones(1,eles)*r_type;