Slightly-less ugly meshplot

function f = meshplot(file, alfa, flag_air, flag_raja)
% Plots a mesh from cim.fedat file or fedat structure

  if nargin == 0
    clc;
    close all;
    file = [fcroot '\ajot\cim.fedat'];
    flag_air = 0;
    flag_raja = 0;
  end,
  if nargin < 2
    alfa = 0;
  end;
  if nargin < 3
    flag_air = 0;
  end;
  if nargin < 4
    flag_raja = 1;
  end;
  flag_air = 1;

  if isstruct(file)
    data = file;
  else
    data = readcimfedat(file,1);
  end;

  % Rotate whole mesh by a given angle
  alfa = alfa*pi/180;
  x = data.cord(:,1);
  y = data.cord(:,2);
  data.cord(:,1) = x*cos(alfa)-y*sin(alfa);
  data.cord(:,2) = y*cos(alfa)+x*sin(alfa);

%   for k = 1 : data.maxel
%     ala = elarea(data.cord(data.nop(k,1:3),1),data.cord(data.nop(k,1:3),2));
%     disp(ala)
%     if ala < 0
%       return
%     end;
%   end;
%   return

%   elx = mask(data.nop(:,1:3), data.cord(:,1))
%   ely = mask(data.nop(:,1:3), data.cord(:,2))
%   cpx = mean(elx, 2);
%   cpy = mean(ely, 2);
%
%   ind = ((cpx.^2 + cpy.^2 > 0.215^2) & (data.imat == 1));
% %
% %   data.nop(ind,:) = [];
%   data.imat(ind,:) = 8;
%
%   iel = [194]
%   data.imat(

%   x1 = data.cord(4000,1);
%   y1 = data.cord(4000,2);
%   x2 = data.cord(3999,1);
%   y2 = data.cord(3999,2);

% imove = 3863
% [kmove xx] = find(data.nop == imove);
% kmove
%
% data.cord(4088,:) = mean(data.cord([5259;4000],:),1);
% data.cord(5260,:) = mean(data.cord([5259;5261],:),1);
% data.cord(5258,:) = mean(data.cord([5257;5259],:),1);
% data.cord(4164,:) = mean(data.cord([4128;5259],:),1);
% data.cord(4059,:) = mean(data.cord([4058;5259],:),1);
%
%
% data.cord(4060,:) = mean(data.cord([3974;3863],:),1);
% data.cord(3864,:) = mean(data.cord([3865;3863],:),1);
% data.cord(4031,:) = mean(data.cord([4032;3863],:),1);
% data.cord(4130,:) = mean(data.cord([4090;3863],:),1);
% data.cord(3862,:) = mean(data.cord([3861;3863],:),1);

%   k = (y2-y1)/(x2-x1);
%
%   data.cord(5259,2) = (y1-k*x1);
%   data.cord(data.ibc(5259),1) = (y1-k*x1);
%
%   data.imat([1907:1912,1942]) = data.imat(1882)
%   data.imat([1883:1888,1913]) = data.imat(1882)

% unique(data.imat)
% return


  % Set rgb colors to each material index
  %setcol = @(ind,c) eval(['col(1+[' num2str(ind) '],:) = ones(' num2str(length(ind)) ',1)*c;']);
  col = zeros(302, 3);
  setcol = @(arg_col, ind, c)( aux_setcol(arg_col, ind, c) );

  col = setcol(col, 0       , [1 1 1]);     % Air
  col = setcol(col, 1       , [1 1 1]*0.6); % Shaft
  col = setcol(col, 2       , [1 1 1]*0.9); % Rotor steel
  col = setcol(col, 3:10    , [1 1 1]*0.7); % Laminations
%   setcol(11      , [0.6 0.6 1]); % Low-permeability layer
  col = setcol(col, 11      , [1 1 0.3]*0.7); % Low-permeability layer
  col = setcol(col, 11      , [1 0 1]); % Low-permeability layer
  col = setcol(col, 32      , [1 1 1]*0.9); % Something
  col = setcol(col, [41 44 61] , [0.4 0.4 1]); % Phase a winding
%   setcol([42 45 62] , [0 1 0]);     % Phase b winding, green
%   setcol([43 46 63] , [1 0 0]);     % Phase c winding, reg
  col = setcol(col, [42 45 62] , [1 0 0]);     % Phase b winding, green
  col = setcol(col, [43 46 63] , [0 1 0]);     % Phase c winding, reg
  col = setcol(col, 47:48   , [1 0.7 0]);   % Rotor windings
  col = setcol(col, 81:92   , [1 0.7 0]);   % Rotor windings
  col = setcol(col, 501:600 , [1 1 0]);     % Rotor bars?
  col = setcol(col, 200:299 , [1 1 0]);     % Rotor bars=
  col = setcol(col, 301     , [1 0 0.8]);   % Permanent magnets
  col = setcol(col, 302     , [1 0 0.5]);   % Permanent magnets
  col = setcol(col, 303     , [1 0 0.2]);   % Permanent magnets

  % 1st order elements
  ne = [3 4 6 8]; % These node numbers are ok
  for in = 1:length(ne)
    if ~any(data.nodel == ne(in)), continue; end;

    % Linear elements
    if ne(in) <= 4
      imats = unique(data.imat);
      if flag_air == 0
        imats(imats == 0) = [];
      end;
      for im = 1:length(imats);
        if flag_raja == 1
          ind = find((data.nodel == ne(in)) & (data.imat == imats(im)));
          p(imats(im)+1) = patch('Vertices', data.cord, 'Faces', data.nop((data.nodel == ne(in)) & (data.imat == imats(im)),1:ne(in)), 'FaceColor', col(imats(im)+1,:)); hold on;
        else
          p(imats(im)+1) = patch('Vertices', data.cord, 'Faces', data.nop((data.nodel == ne(in)) & (data.imat == imats(im)),1:ne(in)), 'FaceColor', col(imats(im)+1,:), 'EdgeColor', 'none'); hold on;
        end;
      end;

    % 2nd order elements, circle arcs needed
    else
      if ne(in) == 6
        ip = [1 4 2; 2 5 3; 3 6 1];
      else
        ip = [1 5 2; 2 6 3; 3 7 4; 4 8 1];
      end;

      kni = find(data.nodel == ne(in));
      for ik = 1: length(kni)
        k = kni(ik);
        if (flag_air == 0) & (data.imat(k) == 0)
          continue;
        end;
        crd = [];
        for i = 1 : size(ip,1);
           x = data.cord(data.nop(k,ip(i,:)),1);
           y = data.cord(data.nop(k,ip(i,:)),2);
           ala = abs(elarea([x, y]));
           if ala > 1e-9
             [xc yc r fi1 fi2] = circlearc(x,y);
             fi = unwrap([fi1 fi2]);
             fi1 = fi(1);
             fi2 = fi(2);
             fii = linspace(fi1,fi2,12); fii([1 end]) = [];
             crd(end+1,:) = data.cord(data.nop(k,i),:);
             for ii = 1 : length(fii)
               crd(end+1,:) = [xc  yc]+r*[cos(fii(ii)) sin(fii(ii))];
             end;
           else
             crd(end+1,:) = data.cord(data.nop(k,i),:);
           end;
        end;
        nop = 1:size(crd,1);
        h = patch('Vertices', crd, 'Faces', nop); hold on;
        set(h, 'FaceColor',col(data.imat(k)+1,:));
        if flag_raja == 0
          set(h, 'EdgeColor','none');
        end;
      end;
    end;
  end;

%   % Plot material limits
  if flag_raja == 0
    sivut = [];
    smat = [];
    for k = 1 : data.maxel
      for i = 1 : 3
        nik1 = data.nop(k,i);
        nik2 = data.nop(k,circind(i+1,3));
        rivi = findrows(sivut, [nik1, nik2]);
        if (k > 1) & length(rivi) == 0
          sivut(end+1,:) = [nik1 nik2];
          smat(end+1) = data.imat(k);
          iraja(length(sivut)) = 1;
        elseif (k == 1)
          sivut(end+1,:) = [nik1 nik2];
          smat(end+1) = data.imat(k);
          iraja(length(sivut)) = 1;
        else
          iraja(rivi) = 0;
          if smat(rivi) ~= data.imat(k);
            iraja(rivi) = 1;
          end;
        end;
      end;
    end;
    iraja = find(iraja)';
    for i = 1 : length(iraja);
      plot(data.cord(sivut(iraja(i),:),1), data.cord(sivut(iraja(i),:),2), 'k-');
    end;
  end;
  set(gcf,'renderer','zbuffer')
  axis equal;

return

%   text(data.cord(:,1), data.cord(:,2), num2strcell(1:data.maxnp))
%   axis equal
%   return

  idir = find(data.ibc == 0);
  iper = find(data.ibc > 0);
  iref = data.ibc(find(data.ibc > 0));
  hold on;
  plot(data.cord(idir, 1), data.cord(idir, 2), 'b.');
  plot(data.cord(iper, 1), data.cord(iper, 2), 'ro');
  plot(data.cord(iref, 1), data.cord(iref, 2), 'go');

  unique(data.cord(iref,1))-unique(data.cord(iper,2))

end

function col = aux_setcol(col, ind, c)

N = size(col, 1);
if max(ind) > N
    col = [col; zeros(max(ind)-N, 3)];
end

col(1+ind, :) = ones( numel(ind) , 1)*c;

end