Index: /issm/trunk-jpl/src/m/classes/planetmesh.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/planetmesh.m	(revision 12913)
+++ /issm/trunk-jpl/src/m/classes/planetmesh.m	(revision 12913)
@@ -0,0 +1,115 @@
+%MESH class definition
+%
+%   Usage:
+%      planetmesh=planetmesh();
+
+classdef planetmesh
+	properties (SetAccess=public) 
+		r                           = NaN;
+		theta                       = NaN;
+		phi                         = NaN
+		elements                    = NaN
+		dimension                   = 0;
+		numberoflayers              = 0;
+		numberofelements            = 0;
+		numberofvertices            = 0;
+		
+		lat                         = NaN
+		long                        = NaN
+
+		vertexconnectivity          = NaN
+		elementconnectivity         = NaN
+		average_vertex_connectivity = 0;
+	end
+	methods
+		function obj = planetmesh(varargin) % {{{
+			switch nargin
+				case 0
+					obj=setdefaultparameters(obj);
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function obj = setdefaultparameters(obj) % {{{
+
+			%the connectivity is the avergaded number of nodes linked to a
+			%given node through an edge. This connectivity is used to initially
+			%allocate memory to the stiffness matrix. A value of 16 seems to
+			%give a good memory/time ration. This value can be checked in
+			%trunk/test/Miscellaneous/runme.m
+			obj.average_vertex_connectivity=25;
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+
+			md = checkfield(md,'planetmesh.r','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
+			md = checkfield(md,'planetmesh.theta','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
+			md = checkfield(md,'planetmesh.phi','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
+			md = checkfield(md,'planetmesh.elements','NaN',1,'>',0,'values',1:md.planetmesh.numberofvertices);
+			if(md.planetmesh.dimension==2),
+				md = checkfield(md,'planetmesh.elements','size',[md.planetmesh.numberofelements 3]);
+			else
+				md = checkfield(md,'planetmesh.elements','size',[md.planetmesh.numberofelements 6]);
+			end
+			if any(~ismember(1:md.planetmesh.numberofvertices,sort(unique(md.planetmesh.elements(:)))));
+				md = checkmessage(md,'orphan nodes have been found. Check the planetmesh outline');
+			end
+			md = checkfield(md,'planetmesh.dimension','values',[2 3]);
+			md = checkfield(md,'planetmesh.numberoflayers','>=',0);
+			md = checkfield(md,'planetmesh.numberofelements','>',0);
+			md = checkfield(md,'planetmesh.numberofvertices','>',0);
+			%no checks for numberofedges lat long and hemisphere
+			if (md.planetmesh.dimension==2),
+				md = checkfield(md,'planetmesh.average_vertex_connectivity','>=',9,'message','''planetmesh.average_vertex_connectivity'' should be at least 9 in 2d');
+			else
+				md = checkfield(md,'planetmesh.average_vertex_connectivity','>=',24,'message','''planetmesh.average_vertex_connectivity'' should be at least 24 in 3d');
+			end
+			md = checkfield(md,'planetmesh.elementconnectivity','size',[md.planetmesh.numberofelements 3],'NaN',1);
+
+			%Solution specific checks
+			switch(solution),
+				case PrognosticSolutionEnum,
+					if md.prognostic.stabilization==3,
+						md = checkfield(md,'planetmesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d planetmeshes');
+					end
+				case TransientSolutionEnum,
+					if md.transient.isprognostic & md.prognostic.stabilization==3,
+						md = checkfield(md,'planetmesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d planetmeshes');
+					end
+				case ThermalSolutionEnum,
+					md = checkfield(md,'planetmesh.dimension','values',3,'message','thermal solution only supported on 3d planetmeshes');
+			end
+		end % }}}
+		function disp(obj) % {{{
+			disp(sprintf('   Mesh:')); 
+
+			disp(sprintf('\n      Elements and vertices:'));
+			fielddisplay(obj,'numberofelements','number of elements');
+			fielddisplay(obj,'numberofvertices','number of vertices');
+			fielddisplay(obj,'elements','index into (x,y,z), coordinates of the vertices');
+			fielddisplay(obj,'r','vertices r coordinate');
+			fielddisplay(obj,'theta','vertices theta coordinate');
+			fielddisplay(obj,'phi','vertices phi coordinate');
+
+			disp(sprintf('\n      Properties:'));
+			fielddisplay(obj,'dimension','planetmesh dimension (2d or 3d)');
+			fielddisplay(obj,'numberoflayers','number of extrusion layers');
+			
+			fielddisplay(obj,'vertexconnectivity','list of vertices connected to vertex_i');
+			fielddisplay(obj,'elementconnectivity','list of vertices connected to element_i');
+			fielddisplay(obj,'average_vertex_connectivity','average number of vertices connected to one vertex');
+
+		end % }}}
+		function marshall(obj,fid) % {{{
+			WriteData(fid,'object',obj,'fieldname','r','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','theta','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','phi','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','elements','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',obj,'fieldname','dimension','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','numberoflayers','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','numberofelements','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','numberofvertices','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','elementconnectivity','format','DoubleMat','mattype',3);
+			WriteData(fid,'object',obj,'fieldname','average_vertex_connectivity','format','Integer');
+		end % }}}
+	end
+end
Index: /issm/trunk-jpl/src/m/model/mesh/planet/mesh_refine_tri4.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/planet/mesh_refine_tri4.m	(revision 12913)
+++ /issm/trunk-jpl/src/m/model/mesh/planet/mesh_refine_tri4.m	(revision 12913)
@@ -0,0 +1,146 @@
+function [ FV ] = mesh_refine_tri4(FV)
+
+% mesh_refine_tri4 - creates 4 triangle from each triangle of a mesh
+%
+% [ FV ] = mesh_refine_tri4( FV )
+%
+% FV.vertices   - mesh vertices (Nx3 matrix)
+% FV.faces      - faces with indices into 3 rows
+%                 of FV.vertices (Mx3 matrix)
+% 
+% For each face, 3 new vertices are created at the 
+% triangle edge midpoints.  Each face is divided into 4
+% faces and returned in FV.
+%
+%        B
+%       /\
+%      /  \
+%    a/____\b       Construct new triangles
+%    /\    /\       [A,a,c]
+%   /  \  /  \      [a,B,b]
+%  /____\/____\     [c,b,C]
+% A	     c	   C    [a,b,c]
+% 
+% It is assumed that the vertices are listed in clockwise order in
+% FV.faces (A,B,C above), as viewed from the outside in a RHS coordinate
+% system.
+% 
+% See also: mesh_refine, sphere_tri, sphere_project
+% 
+
+
+% ---this method is not implemented, but the idea here remains...
+% This can be done until some minimal distance (D) of the mean 
+% distance between vertices of all triangles is achieved.  If
+% no D argument is given, the function refines the mesh once.
+% Alternatively, it could be done until some minimum mean 
+% area of faces is achieved.  As is, it just refines once.
+
+
+% $Revision: 1.1 $ $Date: 2004/11/12 01:32:35 $
+
+% Licence:  GNU GPL, no implied or express warranties
+% History:  05/2002, Darren.Weber_at_radiology.ucsf.edu, created
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+tic;
+fprintf('...refining mesh (tri4)...')
+
+% NOTE
+% The centroid is located one third of the way from each vertex to 
+% the midpoint of the opposite side. Each median divides the triangle 
+% into two equal areas; all the medians together divide it into six 
+% equal parts, and the lines from the median point to the vertices 
+% divide the whole into three equivalent triangles.
+
+% Each input triangle with vertices labelled [A,B,C] as shown
+% below will be turned into four new triangles:
+%
+% Make new midpoints
+% a = (A+B)/2
+% b = (B+C)/2
+% c = (C+A)/2
+%
+%        B
+%       /\
+%      /  \
+%    a/____\b       Construct new triangles
+%    /\    /\       [A,a,c]
+%   /  \  /  \      [a,B,b]
+%  /____\/____\     [c,b,C]
+% A	     c	   C    [a,b,c]
+%
+
+% Initialise a new vertices and faces matrix
+Nvert = size(FV.vertices,1);
+Nface = size(FV.faces,1);
+V2 = zeros(Nface*3,3);
+F2 = zeros(Nface*4,3);
+
+for f = 1:Nface,
+    
+    % Get the triangle vertex indices
+    NA = FV.faces(f,1);
+    NB = FV.faces(f,2);
+    NC = FV.faces(f,3);
+    
+    % Get the triangle vertex coordinates
+    A = FV.vertices(NA,:);
+    B = FV.vertices(NB,:);
+    C = FV.vertices(NC,:);
+    
+    % Now find the midpoints between vertices
+    a = (A + B) ./ 2;
+    b = (B + C) ./ 2;
+    c = (C + A) ./ 2;
+    
+    % Find the length of each median
+    %A2blen = sqrt ( sum( (A - b).^2, 2 ) );
+    %B2clen = sqrt ( sum( (B - c).^2, 2 ) );
+    %C2alen = sqrt ( sum( (C - a).^2, 2 ) );
+    
+    % Store the midpoint vertices, while
+    % checking if midpoint vertex already exists
+    [FV, Na] = mesh_find_vertex(FV,a);
+    [FV, Nb] = mesh_find_vertex(FV,b);
+    [FV, Nc] = mesh_find_vertex(FV,c);
+    
+    % Create new faces with orig vertices plus midpoints
+    F2(f*4-3,:) = [ NA, Na, Nc ];
+    F2(f*4-2,:) = [ Na, NB, Nb ];
+    F2(f*4-1,:) = [ Nc, Nb, NC ];
+    F2(f*4-0,:) = [ Na, Nb, Nc ];
+    
+end
+
+% Replace the faces matrix
+FV.faces = F2;
+
+t=toc; fprintf('done (%5.2f sec)\n',t);
+
+return
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [FV, N] = mesh_find_vertex(FV,vertex)
+
+    Vn = size(FV.vertices,1);
+    Va = repmat(vertex,Vn,1);
+    Vexist = find( FV.vertices(:,1) == Va(:,1) & ...
+                   FV.vertices(:,2) == Va(:,2) & ...
+                   FV.vertices(:,3) == Va(:,3) );
+    if Vexist,
+        if size(Vexist) == [1,1],
+            N = Vexist;
+        else,
+            msg = sprintf('replicated vertices');
+            error(msg);
+        end
+    else
+        FV.vertices(end+1,:) = vertex;
+        N = size(FV.vertices,1);
+    end
+
+return
Index: /issm/trunk-jpl/src/m/model/mesh/planet/planetmesher.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/planet/planetmesher.m	(revision 12913)
+++ /issm/trunk-jpl/src/m/model/mesh/planet/planetmesher.m	(revision 12913)
@@ -0,0 +1,35 @@
+function md=planetmesher(md,varargin)
+%PLANETMESHER - create planet mesh using several packages
+%
+%   This routine creates a planet mesh using several packages: a custom made 
+%   ISSM mesher (called planetmixedmesh.m, method 'mixed') and the planettrimesh
+%   (method 'tria')
+%   where md is a @planet object, varargin is a list of options
+%
+% Usage: md=planetmesh(md,'method','mixed','radius',6378000,'angleresol',1);
+%        md=planetmesh(md,'method','tria','shape','iso','radius',6378000,'refinement',5);
+%
+%
+
+	%recover options
+	options=pairoptions(varargin{:});
+
+	method=getfieldvalue(options,'method','mixed');
+
+	if strcmpi(method,'mixed'),
+		%recover radius and angleresol: 
+		radius=getfieldvalue(options,'radius',6378000); %earth radius as default
+		angleresol=getfieldvalue(options,'angleresol',10); %10 degree resolution
+
+		%call mixed mesh 
+		md=planetmixedmesh(md,radius,angleresol);
+	else
+		%recover radius, shape and level of refinmenet
+		radius=getfieldvalue(options,'radius',6378000); %earth radius as default
+		refinement=getfieldvalue(options,'refinement',5); %refine 5 times
+		shape=getfieldvalue(options,'shape','ico'); 
+
+		%call triangular mesh
+		md=planettrimesh(md,shape,radius,refinement);
+	end
+end
Index: /issm/trunk-jpl/src/m/model/mesh/planet/planetmixedmesh.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/planet/planetmixedmesh.m	(revision 12913)
+++ /issm/trunk-jpl/src/m/model/mesh/planet/planetmixedmesh.m	(revision 12913)
@@ -0,0 +1,99 @@
+function md=planetmixedmesh(md,radius,angleresol)
+%PLANETMIXEDMESH: build 2d shell mesh
+%
+% Usage: md=planetmixedshell(md,radius,angleresol)
+%
+
+conv=pi/180;
+
+r=radius;
+theta=(0:angleresol:360)';
+theta=theta*conv;
+phi=(0:angleresol:180)';
+phi=phi(2:end-1)*conv;
+
+nphi=length(phi);
+nthetha=length(theta);
+
+x=zeros(nphi*nthetha,1);
+y=zeros(nphi*nthetha,1);
+z=zeros(nphi*nthetha,1);
+
+for i=1:nphi,
+	phii=phi(i);
+	x((i-1)*nthetha+1:i*nthetha)=r.*cos(theta).*sin(phii);
+	y((i-1)*nthetha+1:i*nthetha)=r.*sin(theta).*sin(phii);
+	z((i-1)*nthetha+1:i*nthetha)=r.*cos(phii);
+end
+
+quads=zeros((nthetha-1)*(nphi-1),4);
+
+%build elements
+for i=1:nphi-1,
+	for j=1:nthetha-1,
+		count=(i-1)*(nthetha-1)+j;
+	
+		A=(i-1)*nthetha+j;
+		B=(i-1)*nthetha+j+1;
+		C=(i)*nthetha+j+1;
+		D=(i)*nthetha+j;
+		
+		quads(count,:)=[A B C D];
+	end
+end
+
+%now, add polar elements:
+%first north: phi = 0;
+x=[r.*cos(0).*sin(0);x];
+y=[r.*sin(0).*sin(0);y];
+z=[r.*cos(0);z];
+
+%add 1 to quads: 
+quads=quads+1;
+
+%add tria elements:
+trias=zeros(nthetha-1,4);
+
+for i=1:nthetha-1,
+	A=1;
+	B=i+1;
+	C=i+2;
+	trias(i,:)=[A B C NaN];
+end
+
+quads=[trias;quads];
+
+
+% now add south pole: 
+phii=180*conv;
+x=[x;r.*cos(phii).*sin(phii)];
+y=[y;r.*sin(phii).*sin(phii)];
+z=[z;r.*cos(phii)];
+nods=length(x);
+
+%add tria elements:
+trias=zeros(nthetha-1,4);
+
+start=nods-nthetha;
+for i=1:nthetha-1,
+	A=start+i-1;
+	B=start+i;
+	C=nods;
+	trias(i,:)=[A B C NaN];
+end
+
+quads=[quads;trias];
+
+
+md.mesh.elements=quads;
+md.x=x;
+md.y=y;
+md.z=z;
+md.r=sqrt(x.^2+y.^2+z.^2);
+md.theta=acos(z./r);
+md.phi=atan2(y,x);
+
+md.mesh.numberofvertices=length(md.x);
+md.mesh.numberofelements=size(md.mesh.elements,1);
+
+md.mesh.dimension=3;
Index: /issm/trunk-jpl/src/m/model/mesh/planet/planettrimesh.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/planet/planettrimesh.m	(revision 12913)
+++ /issm/trunk-jpl/src/m/model/mesh/planet/planettrimesh.m	(revision 12913)
@@ -0,0 +1,20 @@
+function md=planettrimesh(md,shape,radius,refinement)
+%PLANETTRIMESH: build 2d shell mesh
+%
+% Usage: md=planettrimesh(md,shape,radius,refinement)
+%
+
+results = sphere_tri(shape,refinement,radius);
+md.x=results.vertices(:,1);
+md.y=results.vertices(:,2);
+md.z=results.vertices(:,3);
+md.mesh.elements=results.faces;
+
+md.r=sqrt(md.x.^2+md.y.^2+md.z.^2);
+md.theta=acos(md.z./md.r);
+md.phi=atan2(md.y,md.x);
+
+md.mesh.numberofvertices=length(md.x);
+md.mesh.numberofelements=size(md.mesh.elements,1);
+
+md.mesh.dimension=3;
Index: /issm/trunk-jpl/src/m/model/mesh/planet/runme.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/planet/runme.m	(revision 12913)
+++ /issm/trunk-jpl/src/m/model/mesh/planet/runme.m	(revision 12913)
@@ -0,0 +1,4 @@
+% 5 -> level of refinment
+% 1000 -> radius
+FV = sphere_tri('ico',5,1000); 
+patch('vertices',FV.vertices,'faces',FV.faces,'facecolor',[1 0 0],'edgecolor',[.2 .2 .6]);
Index: /issm/trunk-jpl/src/m/model/mesh/planet/sphere_project.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/planet/sphere_project.m	(revision 12913)
+++ /issm/trunk-jpl/src/m/model/mesh/planet/sphere_project.m	(revision 12913)
@@ -0,0 +1,66 @@
+function V = sphere_project(v,r,c)
+
+% sphere_project - project point X,Y,Z to the surface of sphere radius r
+% 
+% V = sphere_project(v,r,c)
+% 
+% Cartesian inputs:
+% v is the vertex matrix, Nx3 (XYZ)
+% r is the sphere radius, 1x1 (default 1)
+% c is the sphere centroid, 1x3 (default 0,0,0)
+%
+% XYZ are converted to spherical coordinates and their radius is
+% adjusted according to r, from c toward XYZ (defined with theta,phi)
+% 
+% V is returned as Cartesian 3D coordinates
+% 
+
+% $Revision: 1.1 $ $Date: 2004/11/12 01:32:36 $
+
+% Licence:  GNU GPL, no implied or express warranties
+% History:  06/2002, Darren.Weber_at_radiology.ucsf.edu, created
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if ~exist('v','var'),
+    msg = sprintf('SPHERE_PROJECT: No input vertices (X,Y,Z)\n');
+    error(msg);
+end
+
+X = v(:,1);
+Y = v(:,2);
+Z = v(:,3);
+
+if ~exist('c','var'),
+    xo = 0;
+    yo = 0;
+    zo = 0;
+else
+    xo = c(1);
+    yo = c(2);
+    zo = c(3);
+end
+
+if ~exist('r','var'), r = 1; end
+
+% alternate method is to use unit vector of V
+% [ n = 'magnitude(V)'; unitV = V ./ n; ]
+% to change the radius, multiply the unitV
+% by the radius required.  This avoids the
+% use of arctan functions, which have branches.
+
+
+% Convert Cartesian X,Y,Z to spherical (radians)
+theta = atan2( (Y-yo), (X-xo) );
+phi   = atan2( sqrt( (X-xo).^2 + (Y-yo).^2 ), (Z-zo) );
+% do not calc: r = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2);
+
+%   Recalculate X,Y,Z for constant r, given theta & phi.
+R = ones(size(phi)) * r;
+x = R .* sin(phi) .* cos(theta);
+y = R .* sin(phi) .* sin(theta);
+z = R .* cos(phi);
+
+V = [x y z];
+
+return
Index: /issm/trunk-jpl/src/m/model/mesh/planet/sphere_tri.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/planet/sphere_tri.m	(revision 12913)
+++ /issm/trunk-jpl/src/m/model/mesh/planet/sphere_tri.m	(revision 12913)
@@ -0,0 +1,204 @@
+function [FV] = sphere_tri(shape,maxlevel,r,winding)
+
+% sphere_tri - generate a triangle mesh approximating a sphere
+% 
+% Usage: FV = sphere_tri(shape,Nrecurse,r,winding)
+% 
+%   shape is a string, either of the following:
+%   'ico'   starts with icosahedron (most even, default)
+%   'oct'   starts with octahedron
+%   'tetra' starts with tetrahedron (least even)
+%
+%   Nrecurse is int >= 0, setting the recursions (default 0)
+%
+%   r is the radius of the sphere (default 1)
+%
+%   winding is 0 for clockwise, 1 for counterclockwise (default 0).  The
+%   matlab patch command gives outward surface normals for clockwise
+%   order of vertices in the faces (viewed from outside the surface).
+%
+%   FV has fields FV.vertices and FV.faces.  The vertices 
+%   are listed in clockwise order in FV.faces, as viewed 
+%   from the outside in a RHS coordinate system.
+% 
+% The function uses recursive subdivision.  The first
+% approximation is an platonic solid, either an  icosahedron,
+% octahedron or a tetrahedron.  Each level of refinement 
+% subdivides each triangle face by a factor of 4 (see also 
+% mesh_refine).  At each refinement, the vertices are 
+% projected to the sphere surface (see sphere_project).
+% 
+% A recursion level of 3 or 4 is a good sphere surface, if
+% gouraud shading is used for rendering.
+% 
+% The returned struct can be used in the patch command, eg:
+% 
+% % create and plot, vertices: [2562x3] and faces: [5120x3]
+% FV = sphere_tri('ico',4,1);
+% lighting phong; shading interp; figure;
+% patch('vertices',FV.vertices,'faces',FV.faces,...
+%       'facecolor',[1 0 0],'edgecolor',[.2 .2 .6]);
+% axis off; camlight infinite; camproj('perspective');
+% 
+% See also: mesh_refine, sphere_project
+%
+
+
+
+% $Revision: 1.2 $ $Date: 2005/07/20 23:07:03 $
+
+% Licence:  GNU GPL, no implied or express warranties
+% Jon Leech (leech @ cs.unc.edu) 3/24/89
+% icosahedral code added by Jim Buddenhagen (jb1556@daditz.sbc.com) 5/93
+% 06/2002, adapted from c to matlab by Darren.Weber_at_radiology.ucsf.edu
+% 05/2004, reorder of the faces for the 'ico' surface so they are indeed
+% clockwise!  Now the surface normals are directed outward.  Also reset the
+% default recursions to zero, so we can get out just the platonic solids.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+eegversion = '$Revision: 1.2 $';
+fprintf('SPHERE_TRI [v %s]\n',eegversion(11:15)); tic
+
+if ~exist('shape','var') || isempty(shape),
+    shape = 'ico';
+end
+fprintf('...creating sphere tesselation based on %s\n',shape);
+
+% default maximum subdivision level
+if ~exist('maxlevel','var') || isempty(maxlevel) || maxlevel < 0,
+    maxlevel = 0;
+end
+
+% default radius
+if ~exist('r','var') || isempty(r),
+    r = 1;
+end
+
+if ~exist('winding','var') || isempty(winding),
+    winding = 0;
+end
+
+
+% -----------------
+% define the starting shapes
+
+shape = lower(shape);
+
+switch shape,
+case 'tetra',
+    
+    % Vertices of a tetrahedron
+    sqrt_3 = 0.5773502692;
+    
+    tetra.v = [  sqrt_3,  sqrt_3,  sqrt_3 ;   % +X, +Y, +Z  - PPP
+                -sqrt_3, -sqrt_3,  sqrt_3 ;   % -X, -Y, +Z  - MMP
+                -sqrt_3,  sqrt_3, -sqrt_3 ;   % -X, +Y, -Z  - MPM
+                 sqrt_3, -sqrt_3, -sqrt_3 ];  % +X, -Y, -Z  - PMM
+	
+    % Structure describing a tetrahedron
+    tetra.f = [ 1, 2, 3;
+                1, 4, 2;
+                3, 2, 4;
+                4, 1, 3 ];
+    
+    FV.vertices = tetra.v;
+    FV.faces    = tetra.f;
+    
+case 'oct',
+    
+    % Six equidistant points lying on the unit sphere
+    oct.v = [  1,  0,  0 ;  %  X
+              -1,  0,  0 ; 	% -X
+               0,  1,  0 ;  %  Y
+               0, -1,  0 ; 	% -Y
+               0,  0,  1 ; 	%  Z
+               0,  0, -1 ];	% -Z
+	
+    % Join vertices to create a unit octahedron
+    oct.f = [ 1 5 3 ;    %  X  Z  Y  -  First the top half
+              3 5 2 ;    %  Y  Z -X
+              2 5 4 ;    % -X  Z -Y
+              4 5 1 ;    % -Y  Z  X
+              1 3 6 ;    %  X  Y -Z  -  Now the bottom half
+              3 2 6 ;    %  Y  Z -Z
+              2 4 6 ;    % -X  Z -Z
+              4 1 6 ];   % -Y  Z -Z
+    
+    FV.vertices = oct.v;
+    FV.faces    = oct.f;
+    
+case 'ico',
+    
+    % Twelve vertices of icosahedron on unit sphere
+    tau = 0.8506508084; % t=(1+sqrt(5))/2, tau=t/sqrt(1+t^2)
+    one = 0.5257311121; % one=1/sqrt(1+t^2) , unit sphere
+    
+    ico.v( 1,:) = [  tau,  one,    0 ]; % ZA
+    ico.v( 2,:) = [ -tau,  one,    0 ]; % ZB
+    ico.v( 3,:) = [ -tau, -one,    0 ]; % ZC
+    ico.v( 4,:) = [  tau, -one,    0 ]; % ZD
+    ico.v( 5,:) = [  one,   0 ,  tau ]; % YA
+    ico.v( 6,:) = [  one,   0 , -tau ]; % YB
+    ico.v( 7,:) = [ -one,   0 , -tau ]; % YC
+    ico.v( 8,:) = [ -one,   0 ,  tau ]; % YD
+    ico.v( 9,:) = [   0 ,  tau,  one ]; % XA
+    ico.v(10,:) = [   0 , -tau,  one ]; % XB
+    ico.v(11,:) = [   0 , -tau, -one ]; % XC
+    ico.v(12,:) = [   0 ,  tau, -one ]; % XD
+    
+    % Structure for unit icosahedron
+    ico.f = [  5,  8,  9 ;
+               5, 10,  8 ;
+               6, 12,  7 ;
+               6,  7, 11 ;
+               1,  4,  5 ;
+               1,  6,  4 ;
+               3,  2,  8 ;
+               3,  7,  2 ;
+               9, 12,  1 ;
+               9,  2, 12 ;
+              10,  4, 11 ;
+              10, 11,  3 ;
+               9,  1,  5 ;
+              12,  6,  1 ;
+               5,  4, 10 ;
+               6, 11,  4 ;
+               8,  2,  9 ;
+               7, 12,  2 ;
+               8, 10,  3 ;
+               7,  3, 11 ];
+	
+    FV.vertices = ico.v;
+    FV.faces    = ico.f;
+end
+
+
+% -----------------
+% refine the starting shapes with subdivisions
+if maxlevel,
+    
+    % Subdivide each starting triangle (maxlevel) times
+    for level = 1:maxlevel,
+        
+        % Subdivide each triangle and normalize the new points thus
+        % generated to lie on the surface of a sphere radius r.
+        FV = mesh_refine_tri4(FV);
+        FV.vertices = sphere_project(FV.vertices,r);
+        
+        % An alternative might be to define a min distance
+        % between vertices and recurse or use fminsearch
+        
+    end
+end
+
+if winding,
+    fprintf('...returning counterclockwise vertex order (viewed from outside)\n');
+    FV.faces = FV.faces(:,[1 3 2]);
+else
+    fprintf('...returning clockwise vertex order (viewed from outside)\n');
+end
+
+t=toc; fprintf('...done (%6.2f sec)\n\n',t);
+
+return