[src] angefangen Problem in MEX umzusetzen

[src] netzVerfeinerung speichert father2son
[src] funktionen in MEX umgesetzt Auswertung von A fehlt noch

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

diff --git a/src/build_A.cpp b/src/build_A.cpp
new file mode 100644 (file)
index 0000000..021f5b6
--- /dev/null
+++ b/src/build_A.cpp
@@ -0,0 +1,190 @@
+\r
+#include <iostream>\r
+#include <cmath>\r
+#include "mex.h"\r
+#include <stdlib.h>\r
+\r
+#define my_PI 3.141592653589793\r
+\r
+using namespace std;\r
+\r
+int sign(double);\r
+double arsinh(double);\r
+\r
+// sol = g0(p,y,x,l);\r
+double g0(double, double, double, double);\r
+// sol = G00(p,y1,y2,x1,x2,l);\r
+double G00(double, double, double, double, double, double);\r
+// sol = F_par(x1,x2,y1,y2,d1,d2,d3);\r
+double F_par(double, double, double, double, double, double, double);\r
+// sol = F_ort(x1,x2,y1,y2,d1,d2,d3);\r
+//double F_ort(double, double, double, double, double, double, double);\r
+\r
+// sol = quadInt((F_par/F_ort),s1,s2,k1,k2,t1,t2,l1,l2,d1,d2,d3);\r
+double quadInt(double (*)(double,double,double,double,double,double,double),\r
+       double ,double, double ,double ,double ,double ,double ,double ,\r
+       double , double , double );\r
+\r
+void mexFunction(\r
+                int          nlhs,\r
+                mxArray      *plhs[],\r
+                int          nrhs,\r
+                const mxArray *prhs[]\r
+                )\r
+{    \r
+    //sicherheitsabfragen zu Datengroessen\r
+    if (nrhs !=2)\r
+        mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx3))");\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!=3)\r
+        mexErrMsgTxt("expected elements (Mx3)");\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 * 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
+    \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
+    \r
+    double * d = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL));\r
+    \r
+    \r
+    \r
+    //Ausrechnen\r
+for (int j=0;j<em;++j){\r
+    for (int k=0;k<em;++k){\r
+        \r
+     x1[0] = E(j);\r
+        \r
+     // Vorbereiten der DATEN\r
+     //   ej = coordinates(elements(j,:)',:);\r
+     //   ek = coordinates(elements(k,:)',:);\r
+     //   d = ek(1,:) - ej(1,:);\r
+\r
+     // AUSRECHNEN\r
+        A[(k*em)+j] = 1./(4*my_PI);// *\r
+         // quadInt(F_par,\r
+         //   ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));\r
+\r
+    }\r
+ \r
+   \r
+}    \r
+     //Rueckgabe (eventuell zurueck schreiben)     \r
+    \r
+    \r
+ return;   \r
+}\r
+\r
+\r
+       int inline sign(double x)\r
+       {\r
+               return x > 0 ? 1 : (x < 0 ? -1 : 0);\r
+       }\r
+    \r
+    double inline arsinh(double x)\r
+    {\r
+        return log(x+sqrt(x*x+1));\r
+    }\r
+    \r
+    double g0(double p, double y, double x, double l){\r
+               double sol = 0;\r
+\r
+               if(l!=0){\r
+            if(p==0.5)\r
+                sol = (y-x)/2*sqrt((y-x)*(y-x)+\r
+                        l*l)+l*l/(2*arsinh((y-x)/abs(l)));\r
+            else if(p==0)\r
+                sol = y-x;\r
+            else if(p==-0.5)\r
+                sol = arsinh((y-x)/abs(l));\r
+            else if(p==-1)\r
+                sol = atan((y-x)/abs(l));\r
+            else if(p==-1.5)\r
+                sol = (y-x)/((l*l)*sqrt((y-x)*(y-x)+l*l));\r
+            else\r
+                sol = (y-x)*pow((y-x)*(y-x)+l*l,p)+2*p*l*\r
+                        g0(p-1,y,x,l);\r
+               }else{\r
+            if(p==-0.5)\r
+                               sol = sign(y-x)*log(abs(y-x));\r
+                       else\r
+                               sol = (y-x)*pow(abs(y-x),2*p)/(2*p+1);\r
+               }\r
+        return sol;\r
+       }\r
+    \r
+    double G00(double p, double y1,double y2, double x1,double x2, double l) {\r
+               double sol =0;\r
+               if(p==-1.5){\r
+                       if(l==0){\r
+                sol = -sqrt((y1-x1)*(y1-x1)+(y2-x2)*(y2-x2))\r
+                    /((y1-x1)*(y2-x2));\r
+                       }else{\r
+                sol = sign((y1-x1)*(y2-x2))/(2*abs(l))\r
+                    * acos(-2*(y1-x1)*(y1-x1)*(y2-x2)*(y2-x2)\r
+                    /(((y1-x1)*(y1-x1)+l*l)*((y2-x2)*(y2-x2)+l*l))+1);\r
+                       }\r
+\r
+               }else if(p==-0.5){\r
+            if(l!=0)\r
+               sol = 2*p*l*l*G00(p-1,y1,y2,x1,x2,l);  \r
+            if((y1-x1)!=0)\r
+               sol += (y1-x1)*g0(p,y2,x2,sqrt((y1-x1)*(y1-x1)+l*l));\r
+            if((y2-x2)!=0)\r
+               sol += (y2-x2)*g0(p,y1,x1,sqrt((y2-x2)*(y2-x2)+l*l));\r
+            sol /= 2*p+2;\r
+               }else{\r
+            mexErrMsgTxt("no case for p defined");\r
+        }\r
+        \r
+               return sol;\r
+       }\r
+    \r
+    \r
+    double F_par(double x1,double x2,double y1,double y2,double d1,double d2,double d3)\r
+    {\r
+        double sol = (x1-y1-d1)*(x2-y2-d2);\r
+    \r
+        if(sol!=0)\r
+            sol *= G00(-0.5,x1,x2,y1+d1,y2+d2,d3);\r
+\r
+        if((x1-y1-d1)!=0)\r
+            sol -= (x1-y1-d1)*g0(0.5,x1,y1+d1,sqrt((x2-y2-d2)*(x2-y2-d2)+d3*d3));\r
+\r
+        if((x2-y2-d2)!=0)\r
+            sol -= (x2-y2-d2)*g0(0.5,x2,y2+d2,sqrt((x1-y1-d1)*(x1-y1-d1)+d3*d3));\r
+\r
+        double hlp = ((x1-y1-d1)*(x1-y1-d1)+(x2-y2-d2)*(x2-y2-d2)+d3*d3);\r
+        sol +=  1/3.*hlp*sqrt(hlp);\r
+        return sol;\r
+    }\r
+    \r
+    double inline quadInt(double (*f)(double,double,double,double,double,double,double),\r
+            double s1,double s2,double k1,double k2,double t1,double t2,double l1,double l2,\r
+            double d1, double d2, double d3)\r
+    {\r
+              \r
+                  return  f(k1,k2,l1,l2,d1,d2,d3)-f(k1,k2,l1,t2,d1,d2,d3)-f(k1,k2,t1,l2,d1,d2,d3)+f(k1,k2,t1,t2,d1,d2,d3)\r
+    -f(k1,s2,l1,l2,d1,d2,d3)+f(k1,s2,l1,t2,d1,d2,d3)+f(k1,s2,t1,l2,d1,d2,d3)-f(k1,s2,t1,t2,d1,d2,d3)\r
+    -f(s1,k2,l1,l2,d1,d2,d3)+f(s1,k2,l1,t2,d1,d2,d3)+f(s1,k2,t1,l2,d1,d2,d3)-f(s1,l2,t1,t2,d1,d2,d3)\r
+    +f(s1,s2,l1,l2,d1,d2,d3)-f(s1,0,l1,t2,d1,d2,d3)-f(s1,s2,t1,l2,d1,d2,d3)+f(s1,s2,t1,t2,d1,d2,d3);                \r
+     }
\ No newline at end of file

diff --git a/src/g0.m b/src/g0.m
index d576073ea8c9bbcd2f4a51741917fe1edeed8351..be704ebda9f43ed5865786fc1522a4c678150cce 100644 (file)
--- a/src/g0.m
+++ b/src/g0.m
@@ -4,7 +4,8 @@ function sol = g0(p, y, x, lambda)
      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));
+            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)
@@ -14,7 +15,7 @@ function sol = g0(p, y, x, 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);
+            sol = yx*((yx^2 + lambda^2)^p)+2*p*lambda^2*g0(p-1,y,x,lambda);
         end
      else
         if(p==-0.5)

diff --git a/src/refineQuad.m b/src/refineQuad.m
index fffea50f3ab530a51d4469d494524320d0ceceb7..8d31fc9d25829baceaa79a2dab019bb3b2e76e90 100644 (file)
--- a/src/refineQuad.m
+++ b/src/refineQuad.m
@@ -1,7 +1,9 @@
-function [coordinates,elements] = refineQuad(coordinates,elements,type)
+function [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type)
 
 c_loop = size(elements,1);
 
+fa2so = repmat([1:c_loop]',1,4);
+
 for i = 1:c_loop
     c_ele = size(elements,1);
     c_coo = size(coordinates,1);
@@ -27,6 +29,8 @@ 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];
+    
+        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;
@@ -36,6 +40,7 @@ for i = 1:c_loop
         coordinates(c_coo+2,:) = newc;
         
         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;
@@ -45,6 +50,7 @@ for i = 1:c_loop
         coordinates(c_coo+2,:) = newc;
         
         elements(c_ele+1,:) = [c_coo+1,c_coo+2,el(3)];
+        fa2so(i,2)=c_ele+1;
     end
 end
 end

diff --git a/src/test_refine.m b/src/test_refine.m
index 5fd7113aa548970c954f0a63d09c08aa01f4054a..d9cf857351939850a29e38c46a9ed3db26c53acc 100644 (file)
--- a/src/test_refine.m
+++ b/src/test_refine.m
@@ -18,7 +18,7 @@
     end
     
    for i = 1:r_times 
-    [coordinates,elements]=refineQuad(coordinates,elements,marked);
+    [coordinates,elements, f2s]=refineQuad(coordinates,elements,marked);
    end
    
    clear i eles
\ No newline at end of file

diff --git a/src/test_solve1.mat b/src/test_solve1.mat
new file mode 100644 (file)
index 0000000..d220296
Binary files /dev/null and b/src/test_solve1.mat differ