mainbem.m

别人的边界元程序

%% *********** A STATIC POTENTIAL BEM PROGRAM  USING LINEAL ELEMENT************
    %   程序参考. "蔡瑞英等 --边界元法程序设计及工程应用" 1995
    %   采用四节点四边形单元，没有考虑角点影响
%%  variables used in the bem program
global gi ome
global yy
global mtype nord
global x y z
global icod kcd fi tk tq tkt tqt ml
global at bt sum
global xx
global nodbs numbs
global l cx cy cz sol 
global index tp
global b

gi=zeros(10,3);
ome=gi;
yy=zeros(242,1);
nord=zeros(240,4);
x=yy;
y=yy;
z=yy;
icod=nord;
kcd=yy;
fi=yy;
tk=nord;
tq=nord;
tkt=yy;
tqt=yy;
ml=yy;
at=zeros(4,1);
bt=at;
xx=zeros(2928,1);
cx=zeros(20,1);
cy=cx;
cz=cx;
sol=cx;

xk=zeros(4,1);
yk=xk;
zk=xk;
nd=xk;
equ=zeros(242,1);
c=zeros(1,242);
index=zeros(242,1);
b=zeros(240,1);
%%  call input subroutine
xinput;
%% loop acording to each node   ****warning: only one freedom in a node****
for j=1:nodbs
      xp=x(j);
      yp=y(j);
      zp=z(j);
      sum=0;
      yy(j)=0;      %% known varibale in r.h.
      eq(1:nodbs)=0;
      c(1:nodbs)=0;
%%  loop acording to each element
      for k=1:numbs
          ll=0;
          at(:)=0;
          bt(:)=0;
          for i=1:4
              if nord(k,i)==j 
                  ll=i;
              end
              nd(i)=i;                %% what using????
              xk(i)=x(nord(k,i));
              yk(i)=y(nord(k,i));             
              zk(i)=z(nord(k,i));   
          end
          if ll~=0
              impt(ll,nd,xp ,yp,zp,xk,yk,zk);      %% singular integral calculation
          else
              coft(nd,xp,yp,zp,xk,yk,zk);          %% no singular point in the element
          end
%%    form the equation "  AX=B"      
          for i=1:4
              ki=nord(k,i);
              if icod(k,i)==1
                  c(1,ki)=c(1,ki)-bt(i);
                  yy(j)=yy(j)-at(i)*tk(k,i);
              end
              if icod(k,i)==2
                  c(1,ki)=c(1,ki)+at(i);
                  yy(j)=yy(j)+bt(i)*tq(k,i);
              end
              if icod(k,i)==3
                  yy(j)=yy(j)-at(i)*tk(k,i)+bt(i)*tq(k,i);
              end
              sum=sum-at(i);       %% rigid body method to calculate the main coeficiency
          end    
%%  end loop fo each element
      end 
%% 主对角线影响系数的处理
      if mtype==1
      sum=sum+1;    %% if infinte boundary is considered
      end
      if kcd(j)==1 
          c(1,j)=sum;
      else
          yy(j)=yy(j)-sum*fi(j);
      end
%%     using "拟波前法" method to solving the eqations
      equ(1:nodbs)=c(1:nodbs);
      b(j)=yy(j);
      disp('***Now begin solving the equtions, just wating!****')
      xsolver(equ,j,nodbs,1,1.0e-6);
%% end loop of each node
end
%%
yy=b;
%%  将计算得到未知量赋到相应位置上去
for ki=1:numbs
    for kj=1:4
        ij=nord(ki,kj);
        switch icod(ki,kj)
            case 1
                tq(ki,kj)=yy(ij);
            case 2
                tk(ki,kj)=yy(ij);
        end
    end
end
%%   calculate the average potential of boundary nodes
for nj=1:nodbs
    tkt(nj)=0;
    tqt(nj)=0;
    ml(nj)=0;
    for ki=1:numbs
        for kj=1:4
            if nj==nord(ki,kj)
                ml(nj)=ml(nj)+1;
                tkt(nj)=tkt(nj)+tk(ki,kj);
                tqt(nj)=tqt(nj)+tq(ki,kj);
            end
        end
    end
    tkt(nj)=tkt(nj)/ml(nj);
    tqt(nj)=tqt(nj)/ml(nj);
end

%%  calculate the potential of inner nodes
for ip=1:l
    xp=cx(ip);
    yp=cy(ip);
    zp=cx(ip);
    intet(numbs,xp,yp,zp);
    sol(ip)=tp;
end
%%  output the final results
output();