gistfile1.matlab

function [coords, dof, ndof, edof, nelm, ex, ey] = geometry(gridsize)
  if nargin==0
    gridsize = 5; 
  end
  b12 = @(y) .120 - y;
  b16 = @(y) y - .020;
  coords = [];
  for y = 0:gridsize:.070;
    for x = 0:gridsize:.100;
      if(y < .010)
        if(x <= .020 || x >= .080)
          coords = [coords; x y];
        end
      elseif(y < .020)
        if(x <= .040 || x >= .060)
          coords = [coords; x y];
        end
      elseif(y < .060)
        if(x > b16(y) && x < b12(y))
          coords = [coords; x y];
        elseif(x <= b16(y) && b16(y)-x <= gridsize/2)
          coords = [coords; b16(y) y];
        elseif(x >= b12(y) && x-b12(y) <= gridsize/2)
          coords = [coords; b12(y) y];
        end
      elseif(y <= .070)
        if(x >= .040 && x <= .060)
          coords = [coords; x y];
        end
      end      
    end
  end
  ndof=length(coords);
  dof=(1:ndof)';
  % Calculate all triangles possible with delaunay algorithm
  tri = delaunay(coords(:,1),coords(:,2));
  % Go through them all to remove triangles outside the bounds
  remove = [];
  for i=1:length(tri)
    midp = (coords(tri(i,1),:) + coords(tri(i,2),:) + coords(tri(i,3),:))/3;
    % If the midpoint is outside the area, remove the row
    x = midp(1);
    y = midp(2);
    if((y < .010 && (x > .020 && x < .080)) || (y >= .010 && y < .020 && (x > .040 && x < .060)) || (y >= .020 && y < .060 && (x < b16(y) || x > b12(y))) || (y >= .060 && y <= .070 && (x < .040 || x > .060)))
      remove=[remove;i]; % Row to be deleted
    end
  end
  tri(remove,:)=[];
  nelm=length(tri);
  edof = [(1:nelm)' tri];
  [ex,ey]=coordxtr(edof,coords,dof,3);
end