jacob.m

别人的边界元程序

function zxx=jacob(xw,yw,xk,yk,zk)

%% calculate the jacob and normal vector for the problem
global fjacob cosbx cosby cosbz

dx=zeros(4,1);
dy=dx;
%%
dx(1)=-0.25*(1-yw);
dx(2)=0.25*(1-yw);
dx(3)=0.25*(1+yw);
dx(4)=-0.25*(1+yw);

dy(1)=-0.25*(1-xw);
dy(2)=-0.25*(1+xw);
dy(3)=0.25*(1+xw);
dy(4)=0.25*(1-xw);
%%
xd1=dx(1)*xk(1)+dx(2)*xk(2)+dx(3)*xk(3)+dx(4)*xk(4);
yd1=dx(1)*yk(1)+dx(2)*yk(2)+dx(3)*yk(3)+dx(4)*yk(4);
zd1=dx(1)*zk(1)+dx(2)*zk(2)+dx(3)*zk(3)+dx(4)*zk(4);

xd2=dy(1)*xk(1)+dy(2)*xk(2)+dy(3)*xk(3)+dy(4)*xk(4);
yd2=dy(1)*yk(1)+dy(2)*yk(2)+dy(3)*yk(3)+dy(4)*yk(4);
zd2=dy(1)*zk(1)+dy(2)*zk(2)+dy(3)*zk(3)+dy(4)*zk(4);

g1=yd1*zd2-zd1*yd2;
g2=zd1*xd2-xd1*zd2;
g3=xd1*yd2-yd1*xd2;

fjacob=sqrt(g1^2+g2^2+g3^2);
cosbx=g1/fjacob;
cosby=g2/fjacob;
cosbz=g3/fjacob;
%%
zxx=1;
return