meshanellip.m

function [node, face, elem] = meshanellip(c0, rr, tsize, maxvol)
%
% [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)
%

if (nargin < 3)
    error('you must at least provide c0, rr and tsize, see help for details');
end

rr = rr(:)';
if (length(rr) == 1)
    rr = [rr, rr, rr];
elseif (length(rr) == 3)
    % do nothing
elseif (length(rr) ~= 5)
    error('the number of elements for rr is not correct. see help for details');
end

rmax = min(rr(1:3));
if (nargout == 3)
    if (nargin == 3)
        maxvol = tsize * tsize * tsize;
    end
    [node, face, elem] = meshunitsphere(tsize / rmax, maxvol / (rmax * rmax * rmax));
else
    [node, face] = meshunitsphere(tsize / rmax);
end

node = node * diag(rr(1:3));
if (length(rr) == 5)
    theta = rr(4);
    phi = rr(5);
    Rz = [cos(theta) sin(theta) 0; -sin(theta) cos(theta) 0; 0 0 1];
    Rx = [1 0 0; 0 cos(phi) sin(phi); 0 -sin(phi) cos(phi)];
    node = (Rz * Rx * (node'))';
end
node = node + repmat(c0(:)', size(node, 1), 1);