meshanellip

 [node,face,elem]=meshanellip(c0,rr,opt)

 create the surface and tetrahedral mesh of an ellipsoid

 author: Qianqian Fang, <q.fang at neu.edu>

 input: 
   c0:  center coordinates (x0,y0,z0) of the ellipsoid
   rr:  radii of an ellipsoid, 
        if rr is a scalar, this is a sphere with radius rr
        if rr is a 1x3 or 3x1 vector, it specifies the ellipsoid radii [a,b,c]
        if rr is a 1x5 or 5x1 vector, it specifies [a,b,c,theta,phi]
           where theta and phi are the rotation angles along z and x 
           axes, respectively. Rotation is applied before translation.
   tsize: maximum surface triangle size on the sphere
   maxvol: maximu volume of the tetrahedral elements

 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 face element 
   elem: integer array with dimensions of NE x 4, each row represents
         a tetrahedron; if ignored, only produces the surface

 example:
   [node,face,elem]=meshanellip([10,10,-10],[30,20,10,pi/4,pi/4],0.5,0.4);
   plotmesh(node,elem,'x>10');axis equal;

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

0001 function [node,face,elem]=meshanellip(c0,rr,tsize,maxvol)
0002 %
0003 % [node,face,elem]=meshanellip(c0,rr,opt)
0004 %
0005 % create the surface and tetrahedral mesh of an ellipsoid
0006 %
0007 % author: Qianqian Fang, <q.fang at neu.edu>
0008 %
0009 % input:
0010 %   c0:  center coordinates (x0,y0,z0) of the ellipsoid
0011 %   rr:  radii of an ellipsoid,
0012 %        if rr is a scalar, this is a sphere with radius rr
0013 %        if rr is a 1x3 or 3x1 vector, it specifies the ellipsoid radii [a,b,c]
0014 %        if rr is a 1x5 or 5x1 vector, it specifies [a,b,c,theta,phi]
0015 %           where theta and phi are the rotation angles along z and x
0016 %           axes, respectively. Rotation is applied before translation.
0017 %   tsize: maximum surface triangle size on the sphere
0018 %   maxvol: maximu volume of the tetrahedral elements
0019 %
0020 % output:
0021 %   node: node coordinates, 3 columns for x, y and z respectively
0022 %   face: integer array with dimensions of NB x 3, each row represents
0023 %         a surface mesh face element
0024 %   elem: integer array with dimensions of NE x 4, each row represents
0025 %         a tetrahedron; if ignored, only produces the surface
0026 %
0027 % example:
0028 %   [node,face,elem]=meshanellip([10,10,-10],[30,20,10,pi/4,pi/4],0.5,0.4);
0029 %   plotmesh(node,elem,'x>10');axis equal;
0030 %
0031 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0032 %
0033 
0034 if(nargin<3)
0035     error('you must at least provide c0, rr and tsize, see help for details');
0036 end
0037 
0038 rr=rr(:)';
0039 if(length(rr)==1)
0040     rr=[rr,rr,rr];
0041 elseif(length(rr)==3)
0042     % do nothing
0043 elseif(length(rr)~=5)
0044     error('the number of elements for rr is not correct. see help for details');
0045 end
0046 
0047 rmax=min(rr(1:3));
0048 if(nargout==3)
0049     if(nargin==3)
0050     maxvol=tsize*tsize*tsize;
0051     end
0052     [node,face,elem]=meshunitsphere(tsize/rmax,maxvol/(rmax*rmax*rmax));
0053 else
0054     [node,face]=meshunitsphere(tsize/rmax);
0055 end
0056 
0057 node=node*diag(rr(1:3));
0058 if(length(rr)==5)
0059    theta=rr(4);
0060    phi=rr(5);
0061    Rz=[cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0;0 0 1];
0062    Rx=[1 0 0;0 cos(phi) sin(phi);0 -sin(phi) cos(phi)];
0063    node=(Rz*Rx*(node'))';
0064 end
0065 node=node+repmat(c0(:)',size(node,1),1);