Index: /issm/trunk-jpl/src/m/classes/mesh3dsurface.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh3dsurface.py	(revision 20223)
+++ /issm/trunk-jpl/src/m/classes/mesh3dsurface.py	(revision 20223)
@@ -0,0 +1,188 @@
+from MatlabFuncs import *
+from model import *
+from EnumDefinitions import *
+from numpy import *
+
+class mesh3dsurface(object):
+#MESH3DSURFACE class definition
+#
+#   Usage:
+#      mesh3dsurface=mesh3dsurface();
+	def __init__(self,*args): # {{{
+		self.x                           = NaN
+		self.y                           = NaN
+		self.z                           = NaN
+		self.elements                    = NaN
+		self.numberofelements            = 0
+		self.numberofvertices            = 0
+		self.numberofedges               = 0
+
+		self.lat                         = NaN
+		self.long                        = NaN
+		self.r                           = NaN
+
+		self.vertexonboundary            = NaN
+		self.edges                       = NaN
+		self.segments                    = NaN
+		self.segmentmarkers              = NaN
+		self.vertexconnectivity          = NaN
+		self.elementconnectivity         = NaN
+		self.average_vertex_connectivity = 0
+
+		self.extractedvertices           = NaN
+		self.extractedelements           = NaN
+		
+		if len(args)==0:
+			self.setdefaultparameters()
+		elif len(args)==1:
+			self=mesh3dsurface()
+			arg=args[1]
+			fields=fieldnames(arg)
+			for i in range(len(fields)):
+				field=fields[i]
+				if ismember(field,properties('mesh3dsurface')):
+					self.field=arg.field
+		else:
+			raise RuntimeError('constructor not supported')	
+
+		#set defaults
+		self.setdefaultparameters()
+	#}}}
+	def __repr__(self): # {{{
+		string='   2D tria Mesh (horizontal):'
+		string+='\n      Elements and vertices:'
+		string="%s\n%s"%(string,fielddisplay(self,'numberofelements','number of elements'))
+		string="%s\n%s"%(string,fielddisplay(self,'numberofvertices','number of vertices'))
+		string="%s\n%s"%(string,fielddisplay(self,'elements','vertex indices of the mesh elements'))
+		string="%s\n%s"%(string,fielddisplay(self,'x','vertices x coordinate [m]'))
+		string="%s\n%s"%(string,fielddisplay(self,'y','vertices y coordinate [m]'))
+		string="%s\n%s"%(string,fielddisplay(self,'z','vertices z coordinate [m]'))
+		string="%s\n%s"%(string,fielddisplay(self,'lat','vertices latitude [degrees]'))
+		string="%s\n%s"%(string,fielddisplay(self,'long','vertices longitude [degrees]'))
+		string="%s\n%s"%(string,fielddisplay(self,'r','vertices radius [m]'))
+		
+		string="%s\n%s"%(string,fielddisplay(self,'edges','edges of the 2d mesh (vertex1 vertex2 element1 element2)'))
+		string="%s\n%s"%(string,fielddisplay(self,'numberofedges','number of edges of the 2d mesh'))
+
+		string+='\n      Properties:'
+		string="%s\n%s"%(string,fielddisplay(self,'vertexonboundary','vertices on the boundary of the domain flag list'))
+		string="%s\n%s"%(string,fielddisplay(self,'segments','edges on domain boundary (vertex1 vertex2 element)'))
+		string="%s\n%s"%(string,fielddisplay(self,'segmentmarkers','number associated to each segment'))
+		string="%s\n%s"%(string,fielddisplay(self,'vertexconnectivity','list of vertices connected to vertex_i'))
+		string="%s\n%s"%(string,fielddisplay(self,'elementconnectivity','list of vertices connected to element_i'))
+		string="%s\n%s"%(string,fielddisplay(self,'average_vertex_connectivity','average number of vertices connected to one vertex'))
+
+		string+='\n      Extracted model():'
+		string="%s\n%s"%(string,fielddisplay(self,'extractedvertices','vertices extracted from the model()'))
+		string="%s\n%s"%(string,fielddisplay(self,'extractedelements','elements extracted from the model()')) 
+		# }}}
+	def loadobj(self): # {{{
+		# This def is directly called by matlab when a model() selfect is
+		# loaded. Update old properties here
+
+		#2014 Oct. 1st
+		if isstruct(self):
+			oldself=self
+			#Assign property values from struct
+			self=structtoobj(mesh3dsurface(),oldself)
+			if isfield(oldself,'hemisphere'):
+				print ('md.mesh.hemisphere has been automatically converted to EPSG code')
+				if strcmpi(oldself.hemisphere,'n'):
+					self.epsg=3413
+				else:
+					self.epsg=3031
+		# }}}
+	def setdefaultparameters(self): # {{{
+
+		#the connectivity is the averaged 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
+		self.average_vertex_connectivity=25
+		return self
+		# }}}
+	def checkconsistency(self,md,solution,analyses): # {{{
+
+		md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
+		md = checkfield(md,'fieldname','mesh.y','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
+		md = checkfield(md,'fieldname','mesh.z','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
+		md = checkfield(md,'fieldname','mesh.lat','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
+		md = checkfield(md,'fieldname','mesh.long','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
+		md = checkfield(md,'fieldname','mesh.r','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
+		md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',[md.mesh.numberofvertices])
+		md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,3])
+		if any(not ismember([md.mesh.numberofvertices],sort(unique(md.mesh.elements)))):
+			md = checkmessage(md,'orphan nodes have been found. Check the mesh outline')
+		
+		md = checkfield(md,'fieldname','mesh.numberofelements','>',0)
+		md = checkfield(md,'fieldname','mesh.numberofvertices','>',0)
+		md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',9,'message','"mesh.average_vertex_connectivity" should be at least 9 in 2d')
+
+		if (solution==ThermalSolutionEnum()):
+			md = checkmessage(md,'thermal not supported for 2d mesh');
+			
+		return md
+		# }}}
+	def marshall(self,md,fid): # {{{
+		WriteData(fid,'enum',DomainTypeEnum(),'data',StringToEnum(['Domain', domaintype(obj)]),'format','Integer')
+		WriteData(fid,'enum',DomainDimensionEnum(),'data',dimension(obj),'format','Integer')
+		WriteData(fid,'enum',MeshElementtypeEnum(),'data',StringToEnum(elementtype(obj)),'format','Integer')
+		WriteData(fid,'object',obj,'class','mesh','fieldname','x','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',obj,'class','mesh','fieldname','y','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',obj,'class','mesh','fieldname','z','format','DoubleMat','mattype',1)
+		WriteData(fid,'enum',MeshLatEnum(),'data',obj.lat,'format','DoubleMat','mattype',1)
+		WriteData(fid,'enum',MeshLongEnum(),'data',obj.long,'format','DoubleMat','mattype',1)
+		WriteData(fid,'object',obj,'class','mesh','fieldname','r','format','DoubleMat','mattype',1)
+		WriteData(fid,'enum',MeshZEnum(),'data',zeros(obj.numberofvertices,1),'format','DoubleMat','mattype',1)
+		WriteData(fid,'object',obj,'class','mesh','fieldname','elements','format','DoubleMat','mattype',2)
+		WriteData(fid,'object',obj,'class','mesh','fieldname','numberofelements','format','Integer')
+		WriteData(fid,'object',obj,'class','mesh','fieldname','numberofvertices','format','Integer')
+		WriteData(fid,'object',obj,'class','mesh','fieldname','average_vertex_connectivity','format','Integer')
+		WriteData(fid,'object',obj,'class','mesh','fieldname','vertexonboundary','format','DoubleMat','mattype',1)
+		# }}}
+	def domaintype(self): # {{{
+		return '3Dsurface'
+		# }}}
+	def dimension(self): # {{{
+		return 2
+		# }}}
+	def elementtype(self): # {{{
+		return 'Tria'
+		# }}}
+	def processmesh(self,options): # {{{
+	
+		isplanet = 1
+		is2d     = 0
+
+		elements = self.elements
+		x        = self.x
+		y        = self.y
+		z        = self.z
+		return [x, y, z, elements, is2d, isplanet]
+		# }}}
+	def savemodeljs(self,fid,modelname): # {{{
+	
+		fid.write('#s.mesh=new mesh3dsurface()\n'%modelname)
+		writejs1Darray(fid,[modelname, '.mesh.x'],self.x)
+		writejs1Darray(fid,[modelname, '.mesh.y'],self.y)
+		writejs1Darray(fid,[modelname, '.mesh.z'],self.z)
+		writejs2Darray(fid,[modelname, '.mesh.elements'],self.elements)
+		writejsdouble(fid,[modelname, '.mesh.numberofelements'],self.numberofelements)
+		writejsdouble(fid,[modelname, '.mesh.numberofvertices'],self.numberofvertices)
+		writejsdouble(fid,[modelname, '.mesh.numberofedges'],self.numberofedges)
+		writejs1Darray(fid,[modelname, '.mesh.lat'],self.lat)
+		writejs1Darray(fid,[modelname, '.mesh.long'],self.long)
+		writejs1Darray(fid,[modelname, '.mesh.r'],self.r)
+		writejs1Darray(fid,[modelname, '.mesh.vertexonboundary'],self.vertexonboundary)
+		writejs2Darray(fid,[modelname, '.mesh.edges'],self.edges)
+		writejs2Darray(fid,[modelname, '.mesh.segments'],self.segments)
+		writejs2Darray(fid,[modelname, '.mesh.segmentmarkers'],self.segmentmarkers)
+		writejs2Darray(fid,[modelname, '.mesh.vertexconnectivity'],self.vertexconnectivity)
+		writejs2Darray(fid,[modelname, '.mesh.elementconnectivity'],self.elementconnectivity)
+		writejsdouble(fid,[modelname, '.mesh.average_vertex_connectivity'],self.average_vertex_connectivity)
+		writejs1Darray(fid,[modelname, '.mesh.extractedvertices'],self.extractedvertices)
+		writejs1Darray(fid,[modelname, '.mesh.extractedelements'],self.extractedelements)
+
+		# }}}
+	
Index: /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py	(revision 20222)
+++ /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py	(revision 20223)
@@ -4,4 +4,6 @@
 from numpy import *
 from pairoptions import *
+from mesh3dsurface import *
+import subprocess
 
 def gmshplanet(*varargin):
@@ -30,6 +32,6 @@
 
 	#initialize mesh: 
-	#mesh=mesh3dsurface()
-	mesh=mesh2d()
+	mesh=mesh3dsurface()
+	
 
 	#create .geo file:  {{{
@@ -44,42 +46,42 @@
 	fid.write('resolution=%g;\n'%resolution)
 	fid.write('radius=%g;\n'%radius)
-	fid.write('Point(1) = [0.0,0.0,0.0,resolution];\n')
-	fid.write('Point(2) = [radius,0.0,0.0,resolution];\n')
-	fid.write('Point(3) = [0,radius,0.0,resolution];\n')
-	fid.write('Circle(1) = [2,1,3];\n')
-	fid.write('Point(4) = [-radius,0,0.0,resolution];\n')
-	fid.write('Point(5) = [0,-radius,0.0,resolution];\n')
-	fid.write('Circle(2) = [3,1,4];\n')
-	fid.write('Circle(3) = [4,1,5];\n')
-	fid.write('Circle(4) = [5,1,2];\n')
-	fid.write('Point(6) = [0,0,-radius,resolution];\n')
-	fid.write('Point(7) = [0,0,radius,resolution];\n')
-	fid.write('Circle(5) = [3,1,6];\n')
-	fid.write('Circle(6) = [6,1,5];\n')
-	fid.write('Circle(7) = [5,1,7];\n')
-	fid.write('Circle(8) = [7,1,3];\n')
-	fid.write('Circle(9) = [2,1,7];\n')
-	fid.write('Circle(10) = [7,1,4];\n')
-	fid.write('Circle(11) = [4,1,6];\n')
-	fid.write('Circle(12) = [6,1,2];\n')
-	fid.write('Line Loop(13) = [2,8,-10];\n')
-	fid.write('Ruled Surface(14) = [13];\n')
-	fid.write('Line Loop(15) = [10,3,7];\n')
-	fid.write('Ruled Surface(16) = [15];\n')
-	fid.write('Line Loop(17) = [-8,-9,1];\n')
-	fid.write('Ruled Surface(18) = [17];\n')
-	fid.write('Line Loop(19) = [-11,-2,5];\n')
-	fid.write('Ruled Surface(20) = [19];\n')
-	fid.write('Line Loop(21) = [-5,-12,-1];\n')
-	fid.write('Ruled Surface(22) = [21];\n')
-	fid.write('Line Loop(23) = [-3,11,6];\n')
-	fid.write('Ruled Surface(24) = [23];\n')
-	fid.write('Line Loop(25) = [-7,4,9];\n')
-	fid.write('Ruled Surface(26) = [25];\n')
-	fid.write('Line Loop(27) = [-4,12,-6];\n')
-	fid.write('Ruled Surface(28) = [27];\n')
-	fid.write('Surface Loop(29) = [28,26,16,14,20,24,22,18];\n')
-	fid.write('Volume(30) = [29];\n')
-	fid.write('Physical Surface(1) = [28,26,16,14,20,24,22,18];\n')
+	fid.write('Point(1) = {0.0,0.0,0.0,resolution};\n')
+	fid.write('Point(2) = {radius,0.0,0.0,resolution};\n')
+	fid.write('Point(3) = {0,radius,0.0,resolution};\n')
+	fid.write('Circle(1) = {2,1,3};\n')
+	fid.write('Point(4) = {-radius,0,0.0,resolution};\n')
+	fid.write('Point(5) = {0,-radius,0.0,resolution};\n')
+	fid.write('Circle(2) = {3,1,4};\n')
+	fid.write('Circle(3) = {4,1,5};\n')
+	fid.write('Circle(4) = {5,1,2};\n')
+	fid.write('Point(6) = {0,0,-radius,resolution};\n')
+	fid.write('Point(7) = {0,0,radius,resolution};\n')
+	fid.write('Circle(5) = {3,1,6};\n')
+	fid.write('Circle(6) = {6,1,5};\n')
+	fid.write('Circle(7) = {5,1,7};\n')
+	fid.write('Circle(8) = {7,1,3};\n')
+	fid.write('Circle(9) = {2,1,7};\n')
+	fid.write('Circle(10) = {7,1,4};\n')
+	fid.write('Circle(11) = {4,1,6};\n')
+	fid.write('Circle(12) = {6,1,2};\n')
+	fid.write('Line Loop(13) = {2,8,-10};\n')
+	fid.write('Ruled Surface(14) = {13};\n')
+	fid.write('Line Loop(15) = {10,3,7};\n')
+	fid.write('Ruled Surface(16) = {15};\n')
+	fid.write('Line Loop(17) = {-8,-9,1};\n')
+	fid.write('Ruled Surface(18) = {17};\n')
+	fid.write('Line Loop(19) = {-11,-2,5};\n')
+	fid.write('Ruled Surface(20) = {19};\n')
+	fid.write('Line Loop(21) = {-5,-12,-1};\n')
+	fid.write('Ruled Surface(22) = {21};\n')
+	fid.write('Line Loop(23) = {-3,11,6};\n')
+	fid.write('Ruled Surface(24) = {23};\n')
+	fid.write('Line Loop(25) = {-7,4,9};\n')
+	fid.write('Ruled Surface(26) = {25};\n')
+	fid.write('Line Loop(27) = {-4,12,-6};\n')
+	fid.write('Ruled Surface(28) = {27};\n')
+	fid.write('Surface Loop(29) = {28,26,16,14,20,24,22,18};\n')
+	fid.write('Volume(30) = {29};\n')
+	fid.write('Physical Surface(1) = {28,26,16,14,20,24,22,18};\n')
 	fid.write('Physical Volume(2) = 30;\n')
 	fid.close()
@@ -108,8 +110,8 @@
 	#call gmsh
 	if options.exist('refine'):
-		eval('!gmsh -tol 1e-8 -2 sphere.geo -bgm sphere.pos')
+		subprocess.call('gmsh -tol 1e-8 -2 sphere.geo -bgm sphere.pos',shell=True)
 	else:
 		#call gmsh
-		eval('!gmsh -tol 1e-8 -2 sphere.geo')
+		subprocess.call('gmsh -tol 1e-8 -2 sphere.geo',shell=True)
 
 	#import mesh:  {{{
@@ -117,36 +119,42 @@
 
 	#Get Mesh format
-	A=fscanf(fid,'%s',1)
+	A=fid.readline().strip()
 	if not strcmp(A,'$MeshFormat'):
 		raise RuntimeError(['Expecting $MeshFormat (', A, ')'])
 
-	A=fscanf(fid,'%f %i %i',[1, 3])
-	A=fscanf(fid,'%s',1)
+	A=fid.readline().split()
+	A=fid.readline().strip()
 	if not strcmp(A,'$EndMeshFormat'):
 		raise RuntimeError(['Expecting $EndMeshFormat (', A, ')'])
 
 	#Nodes
-	A=fscanf(fid,'%s',1)
+	A=fid.readline().strip()
 	if not strcmp(A,'$Nodes'):
 		raise RuntimeError(['Expecting $Nodes (', A, ')'])
 
-	mesh.numberofvertices=fscanf(fid,'%i',1)
-	A=fscanf(fid,'%i %f %f %f',[4, mesh.numberofvertices])
-	mesh.x = transpose(A[1,:])
-	mesh.y = transpose(A[2,:])
-	mesh.z = transpose(A[3,:])
+	mesh.numberofvertices=int(fid.readline().strip())
+	mesh.x=empty(mesh.numberofvertices)
+	mesh.y=empty(mesh.numberofvertices)
+	mesh.z=empty(mesh.numberofvertices)
+	for i in range(mesh.numberofvertices):
+		A=fid.readline().split()
+		mesh.x[i]=float(A[0])
+		mesh.y[i]=float(A[1])
+		mesh.z[i]=float(A[2])
 
-	A=fscanf(fid,'%s',1)
+	A=fid.readline().strip()
 	if not strcmp(A,'$EndNodes'):
 		raise RuntimeError(['Expecting $EndNodes (', A, ')'])
 
 	#Elements
-	A=fscanf(fid,'%s',1)
+	A=fid.readline().strip()
 	if not strcmp(A,'$Elements'):
 		raise RuntimeError(['Expecting $Elements (', A, ')'])
-	mesh.numberofelements=fscanf(fid,'%i',1)
-	A=fscanf(fid,'%i %i %i %i %i %i %i %i',[8, mesh.numberofelements])
-	mesh.elements=transpose(A[6:8,:]) 
-	A=fscanf(fid,'%s',1)
+	mesh.numberofelements=int(fid.readline().strip())
+	mesh.elements=empty([mesh.numberofelements,3])
+	for i in range(mesh.numberofelements):
+		A=fid.readline().split()
+		mesh.elements[i]=[float(A[5]),float(A[6]),float(A[7])]
+	A=fid.readline().strip()
 	if not strcmp(A,'$EndElements'):
 		raise RuntimeError(['Expecting $EndElements (', A, ')'])
@@ -156,10 +164,9 @@
 	#figure out other fields in mesh3dsurface: 
 	mesh.r=sqrt(mesh.x**2+mesh.y**2+mesh.z**2)
-	mesh.lat = asin(mesh.z/mesh.r)/pi*180
-	mesh.long = atan2(mesh.y,mesh.x)/pi*180
+	mesh.lat = arcsin(mesh.z/mesh.r)/pi*180
+	mesh.long = arctan2(mesh.y,mesh.x)/pi*180
 
 	#erase files: 
-	eval('!rm -rf sphere.geo sphere.msh sphere.pos')
-
+	subprocess.call('rm -rf sphere.geo sphere.msh sphere.pos',shell=True)
 
 	#return mesh: