From: treecity Date: Mon, 18 Apr 2011 16:40:57 +0000 (+0000) Subject: [src] angefangen Problem in MEX umzusetzen X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=8d1d665db03c5e45a8080039673405e89eb5d43e;p=bacc.git [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 index 0000000..021f5b6 --- /dev/null +++ b/src/build_A.cpp @@ -0,0 +1,190 @@ + +#include +#include +#include "mex.h" +#include + +#define my_PI 3.141592653589793 + +using namespace std; + +int sign(double); +double arsinh(double); + +// sol = g0(p,y,x,l); +double g0(double, double, double, double); +// sol = G00(p,y1,y2,x1,x2,l); +double G00(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); +// sol = F_ort(x1,x2,y1,y2,d1,d2,d3); +//double F_ort(double, double, double, double, double, double, double); + +// sol = quadInt((F_par/F_ort),s1,s2,k1,k2,t1,t2,l1,l2,d1,d2,d3); +double quadInt(double (*)(double,double,double,double,double,double,double), + double ,double, double ,double ,double ,double ,double ,double , + double , double , double ); + +void mexFunction( + int nlhs, + mxArray *plhs[], + int nrhs, + const mxArray *prhs[] + ) +{ + //sicherheitsabfragen zu Datengroessen + if (nrhs !=2) + mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx3))"); + if (nlhs >1) + mexErrMsgTxt("has only one output argument"); + + int cm = mxGetM(prhs[0]); + int cn = mxGetN(prhs[0]); + + if(cn!=3) + mexErrMsgTxt("expected coordinates (Nx3)"); + int em = mxGetM(prhs[1]); + int en = mxGetN(prhs[1]); + if(en!=3) + mexErrMsgTxt("expected elements (Mx3)"); + //Vorbereitung der Daten + + plhs[0] = mxCreateDoubleMatrix(em,em,mxREAL); + double * A = mxGetPr(plhs[0]); + double * C = mxGetPr(prhs[0]); + double * E = mxGetPr(prhs[1]); + + double * x1 = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL)); + double * x2 = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL)); + double * x3 = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL)); + double * xn = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL)); + + double * y1 = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL)); + double * y2 = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL)); + double * y3 = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL)); + double * yn = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL)); + + double * d = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL)); + + + + //Ausrechnen +for (int j=0;j 0 ? 1 : (x < 0 ? -1 : 0); + } + + double inline arsinh(double x) + { + return log(x+sqrt(x*x+1)); + } + + double g0(double p, double y, double x, double l){ + double sol = 0; + + if(l!=0){ + if(p==0.5) + sol = (y-x)/2*sqrt((y-x)*(y-x)+ + l*l)+l*l/(2*arsinh((y-x)/abs(l))); + else if(p==0) + sol = y-x; + else if(p==-0.5) + sol = arsinh((y-x)/abs(l)); + else if(p==-1) + sol = atan((y-x)/abs(l)); + else if(p==-1.5) + sol = (y-x)/((l*l)*sqrt((y-x)*(y-x)+l*l)); + else + sol = (y-x)*pow((y-x)*(y-x)+l*l,p)+2*p*l* + g0(p-1,y,x,l); + }else{ + if(p==-0.5) + sol = sign(y-x)*log(abs(y-x)); + else + sol = (y-x)*pow(abs(y-x),2*p)/(2*p+1); + } + return sol; + } + + double G00(double p, double y1,double y2, double x1,double x2, double l) { + double sol =0; + if(p==-1.5){ + if(l==0){ + sol = -sqrt((y1-x1)*(y1-x1)+(y2-x2)*(y2-x2)) + /((y1-x1)*(y2-x2)); + }else{ + sol = sign((y1-x1)*(y2-x2))/(2*abs(l)) + * acos(-2*(y1-x1)*(y1-x1)*(y2-x2)*(y2-x2) + /(((y1-x1)*(y1-x1)+l*l)*((y2-x2)*(y2-x2)+l*l))+1); + } + + }else if(p==-0.5){ + if(l!=0) + sol = 2*p*l*l*G00(p-1,y1,y2,x1,x2,l); + if((y1-x1)!=0) + sol += (y1-x1)*g0(p,y2,x2,sqrt((y1-x1)*(y1-x1)+l*l)); + if((y2-x2)!=0) + sol += (y2-x2)*g0(p,y1,x1,sqrt((y2-x2)*(y2-x2)+l*l)); + sol /= 2*p+2; + }else{ + mexErrMsgTxt("no case for p defined"); + } + + return sol; + } + + + double F_par(double x1,double x2,double y1,double y2,double d1,double d2,double d3) + { + double sol = (x1-y1-d1)*(x2-y2-d2); + + if(sol!=0) + sol *= G00(-0.5,x1,x2,y1+d1,y2+d2,d3); + + if((x1-y1-d1)!=0) + sol -= (x1-y1-d1)*g0(0.5,x1,y1+d1,sqrt((x2-y2-d2)*(x2-y2-d2)+d3*d3)); + + if((x2-y2-d2)!=0) + sol -= (x2-y2-d2)*g0(0.5,x2,y2+d2,sqrt((x1-y1-d1)*(x1-y1-d1)+d3*d3)); + + double hlp = ((x1-y1-d1)*(x1-y1-d1)+(x2-y2-d2)*(x2-y2-d2)+d3*d3); + sol += 1/3.*hlp*sqrt(hlp); + return sol; + } + + double inline quadInt(double (*f)(double,double,double,double,double,double,double), + double s1,double s2,double k1,double k2,double t1,double t2,double l1,double l2, + double d1, double d2, double d3) + { + + 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) + -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) + -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) + +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); + } \ No newline at end of file diff --git a/src/g0.m b/src/g0.m index d576073..be704eb 100644 --- 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 fffea50..8d31fc9 100644 --- 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 5fd7113..d9cf857 100644 --- 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 index 0000000..d220296 Binary files /dev/null and b/src/test_solve1.mat differ