【原创】Gambit->Femlab网格文件

kinggun

偶用Matlab编了一个小文件，可以将Gambit 2.0做出来的2D三角形网格文件转换成Femlab 3.0的网格文件格式，不过可能有点缺陷，请大家指正。[在Polyflow里面可以得到收敛结果，但转换后在Femlab里面不收敛。：（]
  
function mesh = rdmesh(file)
  
if nargin < 1
  error('Function requires one input argument');
elseif ~isstr(file)
  error('Input must be a string representing a filename');
end
  
fid = fopen(file,'r');
if fid==-1
  error('File not found or permission denied');
end
  
% Start of main routine ...
line = fgetl(fid);
while line~=-1
    if ~isempty(findstr(line,'NUMNP'))
        line = fgetl(fid);
        para=sscanf(line,'%f'); % para represents: NUMNP     NELEM     NGRPS    NBSETS     NDFCD     NDFVL
        break;
    end
    line = fgetl(fid);
end
  
while line ~= -1
   % Treat the node coordinates
   if ~isempty(findstr(line,'NODAL COORDINATES'))
      p=[];
      line = fgetl(fid);
      data = sscanf(line,'%f');
      while ~isempty(data)
         p = [p, data(2:3)];
         line = fgetl(fid);
         data = sscanf(line,'%f');
      end
   end
   
   % Treat the triangel element
   if ~isempty(findstr(line,'ELEMENTS/CELLS'))
       t=[];
       line = fgetl(fid);
       data = sscanf(line, '%f');
       while ~isempty(data)
           t = [t, data(4:6)];
           line = fgetl(fid);
           data = sscanf(line, '%f');
       end
       t(4, = 1; % trianlge elements
   end
   
   % Get data of the boundary edge
   if ~isempty(findstr(line,'BOUNDARY CONDITIONS'))
       tmp=[]; % store the repeat boundary edge
       n=1;    % which boundary
       for i=1:para(4)
           line = fgetl(fid);
           line = fgetl(fid);
           data = sscanf(line, '%f');
           while ~isempty(data)
               if data(3)==1
                   tmp=[tmp,[t(1,data(1));t(2,data(1));0;0;n;0;0]];
               elseif data(3)==2
                   tmp=[tmp,[t(2,data(1));t(3,data(1));0;0;n;0;0]];
               else
                   tmp=[tmp,[t(3,data(1));t(1,data(1));0;0;n;0;0]];
               end
               line = fgetl(fid);
               data = sscanf(line, '%f');
           end
           line = fgetl(fid);
           if line == -1
               break;
           end
           n=n+1;
       end
   end
           
   line = fgetl(fid);
end % while
  
fclose(fid);
  
% Treatment of boundary
eps=1E-4;
e=0*tmp;
for ix = 1:para(4)
    n=find(tmp(5,:)==ix);
    s=length(n);
    if abs(p(1,tmp(1,n(1)))-p(1,tmp(2,n(1))))<eps
        tmp(7,n)=1;
        if p(2,tmp(1,n(1)))>p(2,tmp(2,n(1)))
            e(:,n(1):n(s))=fliplr(tmp(:,n(1):n(s)));
            e(1:2,n(1):n(s))=flipud(e(1:2,n(1):n(s)));
        else
            e(:,n(1):n(s))=tmp(:,n(1):n(s));
        end
        len=p(2,e(2,n(s)))-p(2,e(1,n(1)));
        for i=1:s
            e(3,n(i))=(p(2,e(1,n(i)))-p(2,e(1,n(1))))/len;
            e(4,n(i))=(p(2,e(2,n(i)))-p(2,e(1,n(1))))/len;
        end
    else
        tmp(6,n)=1;
        if p(1,tmp(1,n(1)))>p(1,tmp(2,n(1)))
            e(:,n(1):n(s))=fliplr(tmp(:,n(1):n(s)));
            e(1:2,n(1):n(s))=flipud(e(1:2,n(1):n(s)));
        else
            e(:,n(1):n(s))=tmp(:,n(1):n(s));
        end
        len=p(1,e(2,n(s)))-p(1,e(1,n(1)));
        for i=1:s
            e(3,n(i))=(p(1,e(1,n(i)))-p(1,e(1,n(1))))/len;
            e(4,n(i))=(p(1,e(2,n(i)))-p(1,e(1,n(1))))/len;
        end
    end
end
  
% Mesh generate combine with femlab
mesh.p=p; mesh.e=e; mesh.t=t;
mesh.v=[1:para(4); NaN*ones(1,para(4))];
mesh.equiv = zeros(4,0);
save mesh mesh;

simwe

支持原创！