meshacylinder

 [node,face]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
    or
 [node,face,elem]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
 [nplc,fplc]=meshacylinder(c0,c1,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, c1:  cylinder axis end points
   r:   radius of the cylinder; if r contains two elements, it outputs
        a cone trunk, with each r value specifying the radius on each end
   tsize: maximum surface triangle size on the sphere
   maxvol: maximu volume of the tetrahedral elements

         if both tsize and maxvol is set to 0, this function sill return 
         the piecewise-linear-complex (PLC) in the form of the nodes (as node)
         and a cell array (as face).

   ndiv: approximate the cylinder surface into ndiv flat pieces, if 
         ignored, ndiv is set to 20

 output:
   node: node coordinates, 3 columns for x, y and z respectively
   face: integer array with dimensions of NB x 3, each row represents
         a surface mesh triangle
   elem: (optional) integer array with dimensions of NE x 4, each row 
         represents a tetrahedron 

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function [node,face,elem]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
0002 %
0003 % [node,face]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
0004 %    or
0005 % [node,face,elem]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
0006 % [nplc,fplc]=meshacylinder(c0,c1,r,0,0,ndiv);
0007 %
0008 % create the surface and (optionally) tetrahedral mesh of a 3D cylinder
0009 %
0010 % author: Qianqian Fang, <q.fang at neu.edu>
0011 %
0012 % input:
0013 %   c0, c1:  cylinder axis end points
0014 %   r:   radius of the cylinder; if r contains two elements, it outputs
0015 %        a cone trunk, with each r value specifying the radius on each end
0016 %   tsize: maximum surface triangle size on the sphere
0017 %   maxvol: maximu volume of the tetrahedral elements
0018 %
0019 %         if both tsize and maxvol is set to 0, this function sill return
0020 %         the piecewise-linear-complex (PLC) in the form of the nodes (as node)
0021 %         and a cell array (as face).
0022 %
0023 %   ndiv: approximate the cylinder surface into ndiv flat pieces, if
0024 %         ignored, ndiv is set to 20
0025 %
0026 % output:
0027 %   node: node coordinates, 3 columns for x, y and z respectively
0028 %   face: integer array with dimensions of NB x 3, each row represents
0029 %         a surface mesh triangle
0030 %   elem: (optional) integer array with dimensions of NE x 4, each row
0031 %         represents a tetrahedron
0032 %
0033 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0034 %
0035 
0036 if(nargin<3)
0037     error('you must at least provide c0, c1, and r');
0038 end
0039 
0040 if(length(r)==1)
0041     r=[r,r];
0042 end
0043 if(any(r<=0) || all(c0==c1))
0044     error('invalid cylinder parameters');
0045 end
0046 c0=c0(:);
0047 c1=c1(:);
0048 len=sqrt(sum((c0-c1).*(c0-c1)));
0049 % define the axial vector v0 and a perpendicular vector t
0050 v0=c1-c0;
0051 
0052 if(nargin<4)
0053     tsize=min([r,len])/10;
0054 end
0055 
0056 if(nargin<5)
0057     maxvol=tsize^3/5;
0058 end
0059 
0060 % calculate the cylinder end face nodes
0061 if(nargin<6) ndiv=20; end
0062 
0063 dt=2*pi/ndiv;
0064 theta=dt:dt:2*pi;
0065 cx=r(:)*cos(theta);
0066 cy=r(:)*sin(theta);
0067 cx=cx';
0068 cy=cy';
0069 p0=[cx(:,1) cy(:,1) zeros(ndiv,1)];
0070 p1=[cx(:,2) cy(:,2) len*ones(ndiv,1)];
0071 pp=[p0;p1];
0072 no=rotatevec3d(pp,v0)+repmat(c0',size(pp,1),1);
0073 
0074 count=1;
0075 for i=1:ndiv-1
0076    fc{count}={[i i+ndiv i+ndiv+1 i+1],1}; count=count+1;
0077 end
0078 i=ndiv;
0079 fc{count}={[i i+ndiv 1+ndiv 1],1}; count=count+1;
0080 fc{count}={1:ndiv,2};count=count+1;  % bottom inner circle
0081 fc{count}={1+ndiv:2*ndiv,3};count=count+1;  % top inner circle
0082 
0083 if(nargout==2 && tsize==0.0 && maxvol==0.0)
0084     node=no;
0085     face=fc;
0086     return;
0087 end
0088 if(nargin==3)
0089     tsize=len/10;
0090 end
0091 if(nargin<5)
0092     maxvol=tsize*tsize*tsize;
0093 end
0094 [node,elem,face]=surf2mesh(no,fc,min(no),max(no),1,maxvol,[0 0 1],[],0);