coft.m

别人的边界元程序

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

%% calculate the nonsingular element -effection marix
global gi ome
global fjacob cosbx cosby cosbz
global at bt 

at(:)=0;
bt(:)=0;
fn=zeros(4,1);
%%

for ia=1:6
   xl=gi(ia,1);
   wx=ome(ia,1);
   for ja=1:6
       yl=gi(ja,1);
       wy=ome(ja,1);
       fn(1)=0.25*(1-xl)*(1-yl);
       fn(2)=0.25*(1+xl)*(1-yl);
       fn(3)=0.25*(1+xl)*(1+yl);
       fn(4)=0.25*(1-xl)*(1+yl); 
      
       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(xl,yl,xk,yk,zk);
       fl1=wx*wy*fjacob/4/pi;
       fl2=(xqp*cosbx+yqp*cosby+zqp*cosbz)/rq1;
       
       att=-fl1*fl2/rq2;
       btt=fl1/rq1;
       
       for j=1:4
           at(nd(j))=at(nd(j))+fn(j)*att;
           bt(nd(j))=bt(nd(j))+fn(j)*btt;
       end    
%%
   end
end

%%
zx=1;
return