From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Fri, 20 Jan 2012 09:11:41 +0000 (+0000)
Subject: [src] begonnen mit Suche nach eigenartiger Verfeinerung
X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=426cb655e32188eb15ee40a5d60f055088e54ea5;p=bacc.git

[src] begonnen mit Suche nach eigenartiger Verfeinerung
[src] ein paar Kommentare bearbeitet
[src] angefangen refineQuad noch mal zu überarbeiten -> ZIEL: nur verfeinern wenn wirklich nötig

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

diff --git a/src/A_step.m b/src/A_step.m
index 0b20dd4..a76b4c3 100644
--- a/src/A_step.m
+++ b/src/A_step.m
@@ -25,7 +25,7 @@ time = zeros(1,3);
   [coordinates_fine,elements_fine,neigh_fine,f2s]=refineQuad(G_C,G_E,G_N,2);
   time(1) = toc;
 
-  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2))
 
   tic
   data = size(G_E,1);
@@ -50,8 +50,15 @@ time = zeros(1,3);
   end
   time(2) = toc;
   
+  % nur RandElemente Verfeinern
+%   marked = ones(1,size(G_E,1));
+%   marked(find(sum((G_N(:,1:4)==0),2))) = 2;
   
-  marked = mark(x_fine(f2s)',ind,eta,eps);
+  
+  % Adaptiv Verfeinern
+  marked = mark(x_fine(f2s)',ind,eta,eps)
+%   plotShape(coordinates_fine,elements_fine,x_fine);
+%  plotMark(find(marked>1),G_C,G_E);
   
 %   clear 'coordinates_fine' 'elements_fine' 'neigh_fine' 'f2s'
   
diff --git a/src/computeEstSlpMuTilde.m b/src/computeEstSlpMuTilde.m
index 76c6f57..13af47a 100644
--- a/src/computeEstSlpMuTilde.m
+++ b/src/computeEstSlpMuTilde.m
@@ -1,4 +1,10 @@
 function ind=computeEstSlpMuTilde(x_fine,coo,ele,f2s)
+% x_fine - Phi
+% coo - coordinaten
+% ele - elemente
+% f2s - Vater Sohn beziehungen
+
+
 xF2S = x_fine(f2s)';
 
 if(size(xF2S,1)==1)
diff --git a/src/mark.m b/src/mark.m
index 83af46f..01d240e 100644
--- a/src/mark.m
+++ b/src/mark.m
@@ -1,9 +1,10 @@
 function REF = mark(xF2S,ind,theta,eta)
 % function REF = mark(xF2S,ind,theta,eta)
-% xF2S - Father son 
+% xF2S - Father son relation
 % ind - error estimator
 % eta - refine element? (0..1, 1 = All)
 % eps - refine how? (0 = Isotrop)
+% REF - vector with entries [1 : 4] 
 
 if(size(xF2S,1)==1)
   xF2S = xF2S';
@@ -25,18 +26,18 @@ if(theta <1)
 
     ell = find(sum_ind  >= sum_ind(end) * theta,1);
 
-    t1(idx(ell:end)) = 1;
+    t1(idx(ell:end)) = 1;   % Nicht verfeinert
 end
 
 
 %% Wie muss verfeinert werden
-if(eps > 0)
+if(eps > 0) % Horizontal oder Vertikal
     t3 = (eta*abs(Ct(3,:)) >= sqrt(Ct(2,:).^2+Ct(4,:).^2));
     REF(t3) = 3;
     t4 = (eta*abs(Ct(4,:)) >= sqrt(Ct(2,:).^2+Ct(3,:).^2));
     REF(t4) = 4;
 end
-REF(~(t4+t3+t1)) = 2;
+REF(~(t4+t3+t1)) = 2;   % Rest wird Horizontal UND Vertikal geteilt
 
 
 end
diff --git a/src/plotShape.m b/src/plotShape.m
index cf5f42e..d50354d 100644
--- a/src/plotShape.m
+++ b/src/plotShape.m
@@ -72,6 +72,15 @@ if(n)
     end
 end
 
+%     anorm = quadNorm(coordinates,elements);
+%     for idx = 1:eles
+%         current = sum(coordinates(elements(idx,[2,4])',:),1)/2;
+%         current = [current ; current+anorm(idx,:)*ind(idx);coordinates(elements(idx,1)',:)];
+%         plot3(current(:,1),current(:,2),current(:,3),'r'); % Zeichnet Oberflaeche
+%         scatter3(current(2,1),current(2,2),current(2,3),'xr');
+%         hold on
+%     end
+
 xlabel 'x'
 ylabel 'y'
 zlabel 'z'
diff --git a/src/refineQuad.m b/src/refineQuad.m
index 9a2f7d9..4470514 100644
--- a/src/refineQuad.m
+++ b/src/refineQuad.m
@@ -71,12 +71,13 @@ while(1==1)
     break;
   end
   for ele = ref(randperm(length(ref)))  %LOL Zufall bringts
+      
  % # HangingNode Check
     N = G_ref_N(ele,G_ref_N(ele,5:8)==0);
     N2 = N(N~=0); %Nachbarn der Kanten mit nur einem Nachbar
     if(~isempty(N2))
       N3=mod(find((G_ref_N(N2',:)==ele)')-1,4)+5; % ACHTUNG noch mal überprüfen
-      N4=N2(diag(G_ref_N(N2',N3'))~=0);
+      N4=N2(diag(G_ref_N(N2',N3'))~=0)
         if(~isempty(N4))
           for i = N4
             if(i>length(G_ref_t))
@@ -255,10 +256,6 @@ global G_ref_E;
 global G_ref_N;
 global G_ref_f2s;
 
-
-%Innere Beziehungen setzen
-% ...
-
 % plotShape(G_ref_C,G_ref_E);
 
 this = G_ref_N(ele,:);
diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
index a6d212c..79a2dde 100644
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -9,12 +9,9 @@
 #include "slpRectangle.hpp"
 #include "gauss.hpp"
 
-#define quad 2	// Anzahl der Quadratur Auswertungstellen 2^quad
-using namespace std;
-//using namespace slpR;
+#define quad 4	// Anzahl der Quadratur Auswertungstellen 2^quad
 
-//double gauss_w5[] = { 0.1185, 0.2393, 0.2844, 0.2393, 0.1185 };
-//double gauss_n5[] = { 0.0469, 0.2308, 0.5000, 0.7692, 0.9531 };
+using namespace std;
 
 int sign(double);
 //double arsinh(double);
@@ -42,7 +39,7 @@ double inline f_A(double x1, double x2, double y1, double y2, double l) {
 	return 1. / sqrt((x1 - y1) * (x1 - y1) + (x2 - y2) * (x2 - y2) + l * l);
 }
 
-double slpR::g_QY(double p, double y, double x, double l) {
+double inline g_QY(double p, double y, double x, double l) {
 
 	double sol = 0;
 
@@ -57,7 +54,7 @@ double slpR::g_QY(double p, double y, double x, double l) {
 }
 
 //y-x muss != 0 sein
-double slpR::g_AY(double p, double y, double x, double l) {
+double inline g_AY(double p, double y, double x, double l) {
 	//printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);
 
 	double sol = 0;
@@ -88,7 +85,7 @@ double slpR::g_AY(double p, double y, double x, double l) {
 	return sol;
 }
 
-double slpR::G_QY1Y2(double p, double y1, double y2, double x1, double x2,
+double inline G_QY1Y2(double p, double y1, double y2, double x1, double x2,
 		double l) {
 
 	double sol = 0;
@@ -106,7 +103,7 @@ double slpR::G_QY1Y2(double p, double y1, double y2, double x1, double x2,
 	return sol;
 }
 
-double slpR::G_AY2X2(double y1, double y2, double x1, double x2, double d3) {
+double inline G_AY2X2(double y1, double y2, double x1, double x2, double d3) {
 
 //	cout << "("	<< y1 << "," << y2 << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl;
 
@@ -127,7 +124,7 @@ double slpR::G_AY2X2(double y1, double y2, double x1, double x2, double d3) {
 	return sol;
 }
 
-double slpR::G_AY1Y2(double p, double y1, double y2, double x1, double x2,
+double inline G_AY1Y2(double p, double y1, double y2, double x1, double x2,
 		double l) {
 	//	printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);
 	double pt = p + 1.5;
@@ -164,14 +161,14 @@ double slpR::G_AY1Y2(double p, double y1, double y2, double x1, double x2,
 	return sol;
 }
 
-double slpR::Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2,
+double inline Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2,
 		double l) {
 	double sol = 0;
 
 	return sol;
 }
 
-double slpR::F_par(double x1, double x2, double y1, double y2, double d1,
+double inline F_par(double x1, double x2, double y1, double y2, double d1,
 		double d2, double d3) {
 
 	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);
@@ -196,7 +193,7 @@ double slpR::F_par(double x1, double x2, double y1, double y2, double d1,
 	return sol;
 }
 
-double slpR::F_ort(double x1, double x2, double y2, double y3, double d1,
+double inline F_ort(double x1, double x2, double y2, double y3, double d1,
 		double d2, double d3) {
 
 	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y2,y3,d1,d2,d3);
@@ -248,7 +245,7 @@ double slpR::F_ort(double x1, double x2, double y2, double y3, double d1,
 
  }*/
 
-double slpR::apply0Int4(
+double inline apply0Int4(
 		double(*f)(double, double, double, double, double, double, double),
 		double b, double d, double t, double v, double d1, double d2,
 		double d3) {
@@ -264,7 +261,7 @@ double slpR::apply0Int4(
 
 }
 
-double slpR::apply0Int2(
+double inline apply0Int2(
 		double(*f)(double, double, double, double, double, double, double),
 		double b, double d, double t, double v, double d1, double d2,
 		double d3) {
diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp
index fae183d..b7ff935 100644
--- a/src/slpRectangle.hpp
+++ b/src/slpRectangle.hpp
@@ -4,30 +4,30 @@
 namespace slpR {
 
 	// sol = g0(p,y,x,l);
-	double g_AY(double, double, double, double);
+//	double g_AY(double, double, double, double);
 	// sol = g0(p,y,x,l);
-	double g_QY(double, double, double, double);
+//	double g_QY(double, double, double, double);
 
 	// sol = G00(p,y1,y2,x1,x2,l);
-	double G_AY1Y2(double, double, double, double, double, double);
+//	double G_AY1Y2(double, double, double, double, double, double);
 	// sol = G00(p,y1,y2,x1,x2,l);
-	double G_AY2X2(double, double, double, double, double);
+//	double G_AY2X2(double, double, double, double, double);
 	// sol = G00(p,y1,y2,x1,x2,l);
-	double G_QY1Y2(double, double, double, double, double, double);
-	double Gs_AX2Y1Y2(double, double, double, double, double, double);
+//	double G_QY1Y2(double, double, double, double, double, double);
+//	double Gs_AX2Y1Y2(double, double, double, double, double, double);
 
 	// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-	double F_par(double, double, double, double, double, double, double);
+//	double F_par(double, double, double, double, double, double, double);
 	// sol = F_ort(x1,x2,y2,y3,d1,d2,d3);
-	double F_ort(double, double, double, double, double, double, double);
+//	double F_ort(double, double, double, double, double, double, double);
 
 	// sol = quad0Int4((F_par/F_ort),b,d,t,v,d1,d2,d3);
-	double apply0Int4(
-			double(*)(double, double, double, double, double, double, double),
-			double, double, double, double, double, double, double);
-	double apply0Int2(
-			double(*)(double, double, double, double, double, double, double),
-			double, double, double, double, double, double, double);
+//	double apply0Int4(
+//double(*)(double, double, double, double, double, double, double),
+	//double, double, double, double, double, double, double);
+			//	double apply0Int2(
+			//double(*)(double, double, double, double, double, double, double),
+			//double, double, double, double, double, double, double);
 
 	// sol = calcParInt.(b,d,t,v,d1,d2,d3);
 	double calcParIntA(double, double, double, double, double, double, double);
@@ -40,6 +40,8 @@ namespace slpR {
 	double calcOrtIntA(double, double, double, double, double, double, double);
 	double calcOrtIntQX1X2(double, double, double, double, double, double, double);
 
+
+
 	//Voll Analytische Berechnung
 	double calcParIntO0(double, double, double, double, double, double, double,
 			double);