impt.m

别人的边界元程序

function zx=impt(i,nd,xp,yp,zp,xk,yk,zk)

%% calculate the singular element effection matrix
global gi ome
global fjacob cosbx cosby cosbz
global at bt 

fn=zeros(4,1);
xt=zeros(3,1);
yt=zeros(3,1);
fj=xt;
at(:)=0;
bt(:)=0;
%%
m=2;
for ia=1:6
    xl=gi(ia,2);
    wx=ome(ia,2);
    for ja=1:6
        yl=gi(ja,2);
        wy=ome(ja,2);
        
        fl1=tan(0.125*(1+xl)*pi);
        fj(1)=0.125*(1+yl)*pi*(1+fl1^2);
        fj(2)=fj(1);        
        switch i
            case 1
                xt(1)=yl;
                yt(1)=(1+yl)*fl1-1;
                xt(2)=yt(1);
                yt(2)=xt(1);
            case 2
                xt(1)=1-(1+yl)*fl1;
                yt(1)=yl;
                xt(2)=-yt(1);
                yt(2)=-xt(1);               
            case 3
                xt(1)=-yl;
                yt(1)=1-(1+yl)*fl1;
                xt(2)=yt(1);
                yt(2)=xt(1);                
            case 4
                xt(1)=(1+yl)*fl1-1;
                yt(1)=-yl;
                xt(2)=-yt(1);
                yt(2)=-xt(1);
        end
%%  loop acording to subelement--two elment
        for k=1:m
            xw=xt(k);
            yw=yt(k);
            
            fn(1)=0.25*(1-xw)*(1-yw);
            fn(2)=0.25*(1+xw)*(1-yw);
            fn(3)=0.25*(1+xw)*(1+yw);
            fn(4)=0.25*(1-xw)*(1+yw);
       
            xqp=xk(1)*fn(1)+xk(2)*fn(2)+xk(3)*fn(3)+xk(4)*fn(4)-xp;
            yqp=yk(1)*fn(1)+yk(2)*fn(2)+yk(3)*fn(3)+yk(4)*fn(4)-yp;
            zqp=zk(1)*fn(1)+zk(2)*fn(2)+zk(3)*fn(3)+zk(4)*fn(4)-zp;       
            rq2=xqp^2+yqp^2+zqp^2;
            rq1=sqrt(rq2);   
       
             jacob(xw,yw,xk,yk,zk);
             
             fl1=wx*wy*fjacob*fj(k)/4/pi;
             fl2=(xqp*cosbx+yqp*cosby+zqp*cosbz)/rq1;
             
             att=-fl1*fl2/rq2;
             btt=fl1/rq1;
             for j=1:4
                 if j==i 
                 else
                     at(nd(j))=at(nd(j))+fn(j)*att;
                 end
                 bt(nd(j))=bt(nd(j))+fn(j)*btt;
             end
%%
        end
    end
end
%%
zx=1;
return