%   BEM_SIMPLE.m
%   本程序用邊界元方法求解正方形柱體內(nèi)電位分布
%   編程人    沙威(Wei Sha) 安徽大學(xué)(Anhui University) ws108@ahu.edu.cn

clear;clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  1.常數(shù)定義

a=6;          %  正方形長(zhǎng)
N=3;          %  每邊點(diǎn)數(shù)
minstep=a/N;  %  最小離散步長(zhǎng)
TOTAL=N*4;    %  所有點(diǎn)數(shù)
C=1/2;        %  常數(shù)定義
NN=100;       %  積分離散精度 
V_L=300;      %  已知電壓矩陣
xx=a/2;       %  方形內(nèi)部任意一點(diǎn)X坐標(biāo)
yy=a/2;       %  方形內(nèi)部任意一點(diǎn)Y坐標(biāo)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  2.坐標(biāo)定位
%  以方柱左下角為坐標(biāo)原點(diǎn)建立坐標(biāo)系
%  匹配點(diǎn)采用逆時(shí)針?lè)较驈?0,a/6)坐標(biāo)點(diǎn)依次編號(hào)
%  方柱左右兩邊X為常數(shù),方柱上下兩邊Y為常數(shù)

value_b=-minstep/2;   %  下側(cè)初值
value_r=-minstep/2;   %  右側(cè)初值
value_t=a+minstep/2;  %  上測(cè)初值
value_l=a+minstep/2;  %  左側(cè)初值

for i=1:TOTAL;
    
    if (i>0 & i<N+1)          %  下側(cè)
        value_b=value_b+minstep;
        point(1,i)=value_b;
        point(2,i)=0;
    
    elseif (i>N & i<2*N+1)    %  右側(cè)
        value_r=value_r+minstep;
        point(1,i)=a;
        point(2,i)=value_r;
    
    elseif (i>2*N & i<3*N+1)  %  上側(cè)
        value_t=value_t-minstep;
        point(1,i)=value_t;
        point(2,i)=a;

    else                      %  左側(cè)
        value_l=value_l-minstep;
        point(1,i)=0;
        point(2,i)=value_l;
    end;
    
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  3.H矩陣h_st確定

for s=1:TOTAL      %  場(chǎng)點(diǎn)循環(huán)
    for t=1:TOTAL  %  源點(diǎn)循環(huán)  
        
        if (s==t)  %  奇異點(diǎn)處理
            h_st(s,t)=C;
        else
            fieldpoint_x=point(1,s);   %  場(chǎng)點(diǎn)X坐標(biāo)
            currentpoint_x=point(1,t); %  源點(diǎn)X坐標(biāo)
            fieldpoint_y=point(2,s);   %  場(chǎng)點(diǎn)Y坐標(biāo)
            currentpoint_y=point(2,t); %  源點(diǎn)Y坐標(biāo) 
            
            current_x=linspace(currentpoint_x-minstep/2,currentpoint_x+minstep/2,NN); %  X積分變量離散
            current_y=linspace(currentpoint_y-minstep/2,currentpoint_y+minstep/2,NN); %  Y積分變量離散
            
            
            if (t>0 & t<N+1)|(t>2*N & t<3*N+1)  %  上下側(cè)
                
                quad=abs(fieldpoint_y-currentpoint_y)./...
                    ((fieldpoint_x-current_x).^2+(currentpoint_y-fieldpoint_y).^2);
                h_st(s,t)=-(1/(2*pi))*trapz(current_x,quad);
                
                
            else    %  左右側(cè)
                
                quad=abs(fieldpoint_x-currentpoint_x)./...
                    ((fieldpoint_x-currentpoint_x).^2+(current_y-fieldpoint_y).^2);
                h_st(s,t)=-(1/(2*pi))*trapz(current_y,quad);
                
            end;   
        end;
    end;
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  4.K矩陣k_st確定

for s=1:TOTAL       %  場(chǎng)點(diǎn)循環(huán)
    for t=1:TOTAL   %  源點(diǎn)循環(huán)  
        
         if (s==t)  %  奇異點(diǎn)處理
             k_st(s,t)=-(log(minstep/2)-1)*minstep/(2*pi);
             
         else
            fieldpoint_x=point(1,s);   %  場(chǎng)點(diǎn)X坐標(biāo)
            currentpoint_x=point(1,t); %  源點(diǎn)X坐標(biāo)
            fieldpoint_y=point(2,s);   %  場(chǎng)點(diǎn)Y坐標(biāo)
            currentpoint_y=point(2,t); %  源點(diǎn)Y坐標(biāo) 
            
            current_x=linspace(currentpoint_x-minstep/2,currentpoint_x+minstep/2,NN); %  X積分變量離散
            current_y=linspace(currentpoint_y-minstep/2,currentpoint_y+minstep/2,NN); %  Y積分變量離散
            
            
            if ((t>0 & t<N+1)|(t>2*N & t<3*N+1))   %  上下側(cè)
                 
                quad=log( ( (fieldpoint_x-current_x).^2 + ...
                            (currentpoint_y-fieldpoint_y).^2 ).^(1/2) );
                k_st(s,t)=-(1/(2*pi))*trapz(current_x,quad);
                
            else  %  左右側(cè)
                
                quad=log( ( (fieldpoint_x-currentpoint_x).^2 + ...
                            (current_y-fieldpoint_y).^2 ).^(1/2) );
                k_st(s,t)=-(1/(2*pi))*trapz(current_y,quad);
                
            end;
        end;
    end;
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  5.矩陣整序
%   上下側(cè)電荷分布已知
%   左右側(cè)電壓分布已知

H_K1=[h_st(:,[1:N]),-k_st(:,[N+1:2*N]),h_st(:,[2*N+1:3*N]),-k_st(:,[3*N+1:4*N])];
H_K2=[k_st(:,[1:N]),-h_st(:,[N+1:2*N]),k_st(:,[2*N+1:3*N]),-h_st(:,[3*N+1:4*N])];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  6.已知部分電壓和電荷矩陣g_u確定

for u=1:TOTAL;
    
    if ( (u>3*N) & (u<4*N+1) )   %  上側(cè)
        g_u(u)=V_L;    
    else
        g_u(u)=0;
    end;
    
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  7.求解剩下電荷電壓分布并顯示
charge_voltage=(H_K1)^(-1)*H_K2*g_u.';

disp('下側(cè)電位從左到右:');
disp(charge_voltage(1:N));
disp('上側(cè)電位從左到右:')
disp(charge_voltage(3*N:-1:2*N+1));
disp('右側(cè)電荷從上到下:');
disp(charge_voltage(2*N:-1:N+1));
disp('左側(cè)電荷從上到下:');
disp(charge_voltage(3*N+1:4*N));

voltage=[charge_voltage(1:N);g_u(N+1:2*N)';charge_voltage(2*N+1:3*N);g_u(3*N+1:4*N)'];
charge= [g_u(1:N)';charge_voltage(N+1:2*N);g_u(:,[2*N+1:3*N])';charge_voltage(3*N+1:4*N)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  8.方形內(nèi)部測(cè)試點(diǎn)H1矩陣h_st1確定


for t=1:TOTAL  %  源點(diǎn)循環(huán)  
    
    fieldpoint_x=xx;           %  場(chǎng)點(diǎn)X坐標(biāo)
    currentpoint_x=point(1,t); %  源點(diǎn)X坐標(biāo)
    fieldpoint_y=yy;           %  場(chǎng)點(diǎn)Y坐標(biāo)
    currentpoint_y=point(2,t); %  源點(diǎn)Y坐標(biāo) 
    
    current_x=linspace(currentpoint_x-minstep/2,currentpoint_x+minstep/2,NN); %  X積分變量離散
    current_y=linspace(currentpoint_y-minstep/2,currentpoint_y+minstep/2,NN); %  Y積分變量離散
    
    
    if (t>0 & t<N+1)|(t>2*N & t<3*N+1) %  上下側(cè)
        
        quad=abs(fieldpoint_y-currentpoint_y)./...
            ((fieldpoint_x-current_x).^2+(currentpoint_y-fieldpoint_y).^2);
        h_st1(t)=-(1/(2*pi))*trapz(current_x,quad);
        
    else %  左右側(cè)
        
        quad=abs(fieldpoint_x-currentpoint_x)./...
            ((fieldpoint_x-currentpoint_x).^2+(current_y-fieldpoint_y).^2);
        h_st1(t)=-(1/(2*pi))*trapz(current_y,quad);
        
    end;   
end;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  9.方形內(nèi)部測(cè)試點(diǎn)K1矩陣k_st1確定

for t=1:TOTAL  %  源點(diǎn)循環(huán)  
    
    fieldpoint_x=xx;           %  場(chǎng)點(diǎn)X坐標(biāo)
    currentpoint_x=point(1,t); %  源點(diǎn)X坐標(biāo)
    fieldpoint_y=yy;           %  場(chǎng)點(diǎn)Y坐標(biāo)
    currentpoint_y=point(2,t); %  源點(diǎn)Y坐標(biāo) 
    
    current_x=linspace(currentpoint_x-minstep/2,currentpoint_x+minstep/2,NN); %  X積分變量離散
    current_y=linspace(currentpoint_y-minstep/2,currentpoint_y+minstep/2,NN); %  Y積分變量離散
    
    
    if ((t>0 & t<N+1)|(t>2*N & t<3*N+1))   %  上下側(cè)
        
        quad=log( ( (fieldpoint_x-current_x).^2 + ...
            (currentpoint_y-fieldpoint_y).^2 ).^(1/2) );
        k_st1(t)=-(1/(2*pi))*trapz(current_x,quad);
        
    else  %  左右側(cè)
        
        quad=log( ( (fieldpoint_x-currentpoint_x).^2 + ...
            (current_y-fieldpoint_y).^2 ).^(1/2) );
        k_st1(t)=-(1/(2*pi))*trapz(current_y,quad);
        
    end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  10.求解內(nèi)部測(cè)試點(diǎn)電位與解析解
resolve=k_st1*charge-h_st1*voltage;  %  代入離散化泊松公式
show=[xx,yy];

disp('在方柱內(nèi)部電位值');
disp('    x=    y=');
disp(show);
disp('BEM方法為:');
disp(resolve);

analysis=V_L*(a-xx)/a;
disp('解析解為:')
disp(analysis)