meshcylinders.m

function [node,face,elem]=meshcylinders(c0, v, len, varargin)
%
% [node,face]=meshcylinders(c0, v, len, r,tsize,maxvol,ndiv)
%    or
% [node,face,elem]=meshacylinder(c0, v, len, r, tsize,maxvol,ndiv)
% [nplc,fplc]=meshacylinder(c0, v, len,r,0,0,ndiv);
%
% create the surface and (optionally) tetrahedral mesh of a 3D cylinder
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input: 
%   c0, cylinder list axis's starting point
%   v: directional vector of the cylinder
%   len: a scalar or a vector denoting the length of each 
%        cylinder segment along the direction of v
%   tsize, maxvol, ndiv: please see the help for meshacylinder for details
%
% output:
%   node, face, elem: please see the help for meshacylinder for details
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

len=cumsum(len);
[ncyl,fcyl]=meshacylinder(c0,c0+v*len(1),varargin{:});

for i=2:length(len)
   [ncyl1,fcyl1]=meshacylinder(c0+v*len(i-1),c0+v*len(i),varargin{:});
   fcyl1=cellfun(@(x) {x{1}+size(ncyl,1),x{2}}, fcyl1, 'UniformOutput', false);
   ncyl=[ncyl; ncyl1];
   if(i==1)
       fcyl1=fcyl1(1:end-1);
   else
       fcyl1={fcyl1{1:end-2},fcyl1{end}};
   end
   fcyl={fcyl{:}, fcyl1{:}};
end

[ncyl,I,J]=unique(round(ncyl*1e10),'rows');
ncyl=ncyl*1e-10;
fcyl=cellfun(@(x) {J(x{1})',x{2}}, fcyl, 'UniformOutput', false);

tsize=varargin{2};
maxvol=varargin{3};

if(nargout==2 && tsize==0.0 && maxvol==0.0)
    node=ncyl;
    face=fcyl;
    return;
end
if(nargin==3)
    tsize=len/10;
end
if(nargin<5)
    maxvol=tsize*tsize*tsize;
end

centroid=cumsum([0 len(1:end-1)])+len/2;   % define the centroids of each cylinder segment
seeds=repmat(c0(:)',length(len),1)+repmat(v(:)',length(len),1).*repmat(centroid(:),1,3);
[node,elem,face]=surf2mesh(no,fc,min(no),max(no),1,maxvol,seeds,[],0);