potential2d.m

本程序为边界元的原程序

function potential2D( file_in )
global G H dfi du1 du2 dq1 dq2 pot
%this programme solves 2D potential problem using constant boundary element
%      potential2D( filename )
%  输入参数： 
%      file_in  ---------- 边界元模型文件
if nargin < 1
        file_in = 'potential2D.dat' ;
    end
    
    % 检查文件是否存在
    if exist( file_in ) == 0
        disp( sprintf( '错误：文件 %s 不存在', file_in ) )
        disp( sprintf( '程序终止' ) )
        return ;
    end
     % 根据输入文件名称生成输出文件名称
    [path_str,name_str,ext_str] = fileparts( file_in ) ;
    ext_str_out = '.mat' ;
    file_out = fullfile( path_str, [name_str, ext_str_out] ) ;
    
    % 检查输出文件是否存在
    if exist( file_out ) ~= 0 
        answer = input( sprintf( '文件 %s 已经存在，是否覆盖? ( Yes / [No] ):  ', file_out ), 's' ) ;
        if length( answer ) == 0
            answer = 'no' ;
        end
        
        answer = lower( answer ) ;
        if answer ~= 'y' | answer ~= 'yes' 
            disp( sprintf( '请使用另外的文件名，或备份已有的文件' ) ) ;
            disp( sprintf( '程序终止' ) );
            return ; 
        end
    end
    
    BemModel( file_in ) ;       % 定义有限元模型
    SolveModel;                 %求解边界元模型
    SolveFunction(G,dfi);       %求解AX=F
    SolveInternal;              %求解内点的位势及位势梯度
    DisplayResults ;            % 把计算结果显示出来 
return

function BemModel( file_in)
%  定义有限元模型
%  输入参数：
%      file_in ---- 有限模型数据文件
%  返回值：
%      无

%  全局变量及名称
%cx:内点的x1坐标
%cy:内点的x2坐标
%x:边界元端点x1坐标值的一维数组
%y:边界元端点x2坐标值的一维数组
%kode:指定节点已知边界条件类型的一维数组，对于J边界节点上规定kode（j）＝0则位势已知，kode（j）＝1则位势梯度已知
%fi:已知边界条件的一维数组，j节点的已知边界条件进入于fi（j）
global cx cy x y kode fi N L 
% 打开文件
    fid = fopen( file_in, 'r' ) ;
    
    %边界点和内点的坐标
    N= fscanf( fid, '%d', 1 ) ;  %边界点数目
    L= fscanf( fid, '%d', 1 ) ;  %内点数目
    x = zeros( N, 1 ) ;
    y = zeros( N, 1 ) ;
    cx = zeros( L, 1 ) ;
    cy = zeros( L, 1 ) ;
    for i = 1:N                  %边界点坐标
        x(i) = fscanf( fid, '%f', 1 ) ;
        y(i)=fscanf(fid, '%f',1);
    end
    for i = 1:L                  %内点坐标
        cx(i) = fscanf( fid, '%f', 1 ) ;
        cy(i)=fscanf(fid, '%f',1);
    end
    
   %各边界点的边界条件（位势及位势导数的值）
    kode = zeros( N, 1 ) ;
    fi = zeros( N, 1 ) ;
   for j=1:N
       kode(j)=fscanf(fid,'%d',1);
       fi(j)=fscanf(fid,'%f',1);
   end
  return
  
  function SolveModel   %求解边界元模型
  %xm:   节点的x1坐标(单元中点x1坐标)
  %ym:   节点的x2坐标(单元中点x2坐标)
  %G、H: HU=GQ
  %dfi: 等号右侧的一维数组，求解方程后存储未知的u和q

  global x y xm ym G H du1 du2 dq1 dq2 fi dfi kode NX N L
 
  %计算节点坐标并存入xm、ym中
  x(N+1)=x(1);    
  y(N+1)=y(1);
  for i=1:N
      xm(i)=(x(i)+x(i+1))/2;
      ym(i)=(y(i)+y(i+1))/2;
  end
  
  %生成和装配G、H矩阵
  for i=1:N
      for j=1:N
          KK=j+1;
          if i-j~=0
             NonDiag(i,j,xm(i),ym(i),x(j),y(j),x(KK),y(KK),0);%求非对角线的G(i,j),H(i,j).这里求得的G，H都没有×2pi。i,j自己加的，H（i,j）,G(i,j),dq1,dq2,du1,du2没加
          end
          if i-j==0
              Duijiaoxian(i,j,x(j),y(j),x(KK),y(KK));%求对角线的G(i,j)，这里的G是实际的×2pi
              H(i,j)=pi;                             %对角线上的H实际为1/2，这里×2pi
          end
      end
  end
    
  
  %再通过已知边界条件移项形成A、F(AX=F)A存储于G，F存于dfi
  for j=1:N
      if kode(j)==1                           %当位势梯度已知，即位势未知，将未知边值及其影响系数放于左侧，即放入G中
          for i=1:N
              CH=G(i,j);
              G(i,j)=-H(i,j);
              H(i,j)=-CH;
          end
      end
  end
  dfi=H*fi;
    
  return
 
  
 
     
    function SolveFunction(A,B) %高斯消去法，得到的B即为所需求的解，A并非为对角阵
  
    global   NX N G dfi
    TOL=1e-6;
    n1=N-1;
    for k=1:n1
        k1=k+1;
        c=A(k,k);
  
        if abs((abs(c)-TOL))<=TOL  %如果对角线为0需换行
            for j=k1:N
                if (abs(A(j,k)-TOL))>0
                    for l=k:N
                        c=A(k,l);
                        A(k,l)=A(j,l);
                        A(j,l)=c;
                    end
                    c=B(k);
                    B(k)=B(j);
                    B(j)=c;
                    c=A(k,k);
                   break
                end
            end
            if (abs(A(j,k)-TOL))<=0
                fprintf('singularity in row\n');
                return
                %加个终止整个程序的命令
            end
        end
        %divide row by diagonal coefficient
   
        c=A(k,k);
        for j=k1:N
            A(k,j)=A(k,j)/c;
        end
            B(k)=B(k)/c;
        %eliminate unknown x(k) from row I
        for i=k1:N
            c=A(i,k);
            for j=k1:N
                A(i,j)=A(i,j)-c*A(k,j);
            end
            B(i)=B(i)-c*B(k);
        end
    end
       
        %compute last unknown
        if (abs(A(N,N))-TOL)~=0
            B(N)=B(N)/A(N,N);
        end
        if (abs(A(N,N))-TOL)==0
            fprintf('singularity in row\n');
                %加个终止整个程序的命令
        end
        
        %反算所有未知量
        for i=1:n1
            k=N-i;
            k1=k+1;
            for j=k1:N
                B(k)=B(k)-A(k,j)*B(j);
            end
        end
        
        dfi=B;
        return
        
        function SolveInternal  %求内点
        global cx cy x y kode fi N L G H dfi du1 du2 dq1 dq2 flux1 flux2 pot
        
        %重新排列fi、dfi，使得fi中存入位势而dfi中存入位势梯度
        for i=1:N
            if kode(i)==1
                ch=fi(i);
                fi(i)=dfi(i);
                dfi(i)=ch;
            end
        end
        
        if L==0
            return
            fprintf('无内点'); 
        end
        for k=1:L     %k为内点
            pot(k)=0;
            flux1(k)=0;
            flux2(k)=0;
            for j=1:N  %内点到各个边界点
                kk=j+1;
                NonDiag(k,j,cx(k),cy(k),x(j),y(j),x(kk),y(kk),1);
                pot(k)=pot(k)+dfi(j)*G(k,j)-fi(j)*H(k,j);
                flux1(k)=flux1(k)+dfi(j)*du1-fi(j)*dq1;
                flux2(k)=flux2(k)+dfi(j)*du2-fi(j)*dq2;
            end
            pot(k)=pot(k)/(2*pi);
            flux1(k)=flux1(k)/(2*pi);
            flux2(k)=flux2(k)/(2*pi);
        end
      
        return
       
  
  function DisplayResults
%  显示计算结果
%  输入参数：
%     无
%  返回值：
%     无
global x y N L kode fi dfi pot flux1 flux2 cx cy xm ym

 fprintf( '\n\n\n' ) ;
    fprintf( '                        表1 边界点位势及位势梯度\n' ) ;
    fprintf( '================================================================================================\n' ) ;
    fprintf( '    坐标X                 坐标Y                      位势                  位势梯度 \n' ) ;
    fprintf( '------------------------------------------------------------------------------------------------\n' ) ;
    for i=1:N
        fprintf( '%.4f                  %.4f                      %.2f                  %.2f \n',xm(i), ym(i),fi(i),dfi(i) ) ;    
    end
    fprintf( '================================================================================================\n\n' ) ;
    
   fprintf( '\n\n\n' ) ;
    fprintf( '                        表2 内点位势及位势梯度\n' ) ;
    fprintf( '============================================================================================================\n' ) ;
    fprintf( '    坐标X              坐标Y                 位势                位势梯度1                位势梯度2 \n' )  ;
    fprintf( '------------------------------------------------------------------------------------------------------------\n' ) ;
    for i=1:L
        fprintf( '%.4f                %.4f                 %.2f               %.3f                   %.5f \n',cx(i),cy(i),pot(i),flux1(i),flux2(i) ) ;    
    end
    fprintf( '=============================================================================================================\n\n' ) 
  return