Index: issm/trunk/src/m/mesh/bamg.py
===================================================================
--- issm/trunk/src/m/mesh/bamg.py	(revision 21341)
+++ issm/trunk/src/m/mesh/bamg.py	(revision 22758)
@@ -1,5 +1,7 @@
 import os.path
 import numpy as  np
-from mesh2d import mesh2d
+from mesh2d import *
+from mesh2dvertical import *
+from mesh3dsurface import *
 from collections import OrderedDict
 from pairoptions import pairoptions
@@ -20,4 +22,8 @@
 
 	   - domain :            followed by an ARGUS file that prescribes the domain outline
+	   - holes :             followed by an ARGUS file that prescribes the holes
+     - subdomains :        followed by an ARGUS file that prescribes the list of
+	                         subdomains (that need to be inside domain)
+
 	   - hmin :              minimum edge length (default is 10^-100)
 	   - hmax :              maximum edge length (default is 10^100)
@@ -37,5 +43,4 @@
 	                         1 -> use Green formula
 	   - KeepVertices :      try to keep initial vertices when adaptation is done on an existing mesh (default 1)
-	   - MaxCornerAngle :    maximum angle of corners in degree (default is 10)
 	   - maxnbv :            maximum number of vertices used to allocate memory (default is 10^6)
 	   - maxsubdiv :         maximum subdivision of exisiting elements (default is 10)
@@ -49,5 +54,4 @@
 	   - power :             power applied to the metric (default is 1)
 	   - splitcorners :      split triangles whuch have 3 vertices on the outline (default is 1)
-	   - geometricalmetric : take the geometry into account to generate the metric (default is 0)
 	   - verbose :           level of verbosity (default is 1)
 
@@ -74,26 +78,40 @@
 	bamg_mesh=bamgmesh()
 
+	subdomain_ref = 1
+	hole_ref = 1
+
 	# Bamg Geometry parameters {{{
 	if options.exist('domain'):
-
 		#Check that file exists
 		domainfile=options.getfieldvalue('domain')
-		if not os.path.exists(domainfile):
-			raise IOError("bamg error message: file '%s' not found" % domainfile)
-		domain=expread(domainfile)
+		if type(domainfile) == str:
+			if not os.path.exists(domainfile):
+				raise IOError("bamg error message: file '%s' not found" % domainfile)
+			domain=expread(domainfile)
+		else:
+			domain=domainfile
 
 		#Build geometry 
 		count=0
 		for i,domaini in enumerate(domain):
-
 			#Check that the domain is closed
 			if (domaini['x'][0]!=domaini['x'][-1] or domaini['y'][0]!=domaini['y'][-1]):
 				raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed")
 
-			#Checks that all holes are INSIDE the principle domain outline
+			#Checks that all holes are INSIDE the principle domain outline princial domain should be index 0
 			if i:
 				flags=ContourToNodes(domaini['x'],domaini['y'],domainfile,0)[0]
 				if np.any(np.logical_not(flags)):
 					raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
+
+			#Check orientation
+			nods=domaini['nods']-1#the domain are closed 1=end;
+
+			test = np.sum((domaini['x'][1:nods+1] - domaini['x'][0:nods])*(domaini['y'][1:nods+1] + domaini['y'][0:nods]))
+			if (i==0 and test>0) or (i>0 and test<0):
+				print('At least one contour was not correctly oriented and has been re-oriented')
+				domaini['x'] = np.flipud(domaini['x'])
+				domaini['y'] = np.flipud(domaini['y'])
+
 
 			#Add all points to bamg_geometry
@@ -102,13 +120,91 @@
 			bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
 			if i:
-				bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,1]))
-
-			# bamg_geometry.Vertices=np.hstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods].reshape(-1),domaini['y'][0:nods].reshape(-1),np.ones((nods))))))
-			# bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.hstack((np.arange(count+1,count+nods+1).reshape(-1,),np.hstack((np.arange(count+2,count+nods+1),count+1)).reshape(-1,),1.*np.ones((nods,))))))
-			# if i:
-			# 	bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,1]))
-
+				bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,-subdomain_ref]))
+				subdomain_ref = subdomain_ref+1;
+			else:
+				bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,0]))
+			
 			#update counter
 			count+=nods
+
+		#Deal with domain holes
+		if options.exist('holes'):
+			holesfile=options.getfieldvalue('holes')
+			if type(holesfile) == str:
+				if not os.path.exists(holesfile):
+					raise IOError("bamg error message: file '%s' not found" % holesfile)
+				holes=expread(holesfile)
+			else:
+				holes=holesfile
+
+			#Build geometry 
+			for i,holei in enumerate(holes):
+				#Check that the hole is closed
+				if (holei['x'][0]!=holei['x'][-1] or holei['y'][0]!=holei['y'][-1]):
+					raise RuntimeError("bamg error message: all contours provided in ''hole'' should be closed")
+
+				#Checks that all holes are INSIDE the principle domain outline princial domain should be index 0
+				flags=ContourToNodes(holei['x'],holei['y'],domainfile,0)[0]
+				if np.any(np.logical_not(flags)):
+					raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
+
+				#Check orientation
+				nods=holei['nods']-1#the hole are closed 1=end;
+				test = np.sum((holei['x'][1:nods+1] - holei['x'][0:nods])*(holei['y'][1:nods+1] + holei['y'][0:nods]))
+				if test<0:
+					print('At least one hole was not correctly oriented and has been re-oriented')
+					holei['x'] = np.flipud(holei['x'])
+					holei['y'] = np.flipud(holei['y'])
+
+				#Add all points to bamg_geometry
+				nods=holei['nods']-1    #the hole are closed 0=end
+				bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((holei['x'][0:nods],holei['y'][0:nods],np.ones((nods)))).T))
+				bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
+				#update counter
+				count+=nods
+
+		#And subdomains
+		if options.exist('subdomains'):
+			subdomainfile=options.getfieldvalue('subdomains')
+			if type(subdomainfile) == str:
+				if not os.path.exists(subdomainfile):
+					raise IOError("bamg error message: file '%s' not found" % subdomainfile)
+				subdomains=expread(subdomainfile)
+			else:
+				subdomains=subdomainfile
+
+			#Build geometry 
+			for i,subdomaini in enumerate(subdomains):
+				#Check that the subdomain is closed
+				if (subdomaini['x'][0]!=subdomaini['x'][-1] or subdomaini['y'][0]!=subdomaini['y'][-1]):
+					raise RuntimeError("bamg error message: all contours provided in ''subdomain'' should be closed")
+
+				#Checks that all subdomains are INSIDE the principle subdomain outline princial domain should be index 0
+				if i:
+					flags=ContourToNodes(subdomaini['x'],subdomaini['y'],domainfile,0)[0]
+					if np.any(np.logical_not(flags)):
+						raise RuntimeError("bamg error message: All subdomains should be strictly inside the principal subdomain")
+
+				#Check orientation
+				nods=subdomaini['nods']-1#the subdomain are closed 1=end;
+
+				test = np.sum((subdomaini['x'][1:nods+1] - subdomaini['x'][0:nods])*(subdomaini['y'][1:nods+1] + subdomaini['y'][0:nods]))
+				if test>0:
+					print('At least one subcontour was not correctly oriented and has been re-oriented')
+					subdomaini['x'] = np.flipud(subdomaini['x'])
+					subdomaini['y'] = np.flipud(subdomaini['y'])
+
+				#Add all points to bamg_geometry
+				nods=subdomaini['nods']-1    #the subdomain are closed 0=end
+				bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((subdomaini['x'][0:nods],subdomaini['y'][0:nods],np.ones((nods)))).T))
+				bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
+				#update counter
+				count+=nods
+			
+		if options.getfieldvalue('vertical',0):
+			if np.size(options.getfieldvalue('Markers',[]))!=np.size(bamg_geometry.Edges,0):
+				raise RuntimeError('for 2d vertical mesh, ''Markers'' option is required, and should be of size ' + str(np.size(bamg_geometry.Edges,0)))
+			if np.size(options.getfieldvalue('Markers',[]))==np.size(bamg_geometry.Edges,0):
+				bamg_geometry.Edges[:,2]=options.getfieldvalue('Markers')
 
 		#take care of rifts
@@ -303,5 +399,4 @@
 	bamg_options['hVertices']=options.getfieldvalue('hVertices',np.empty((0,1)))
 	bamg_options['KeepVertices']=options.getfieldvalue('KeepVertices',1)
-	bamg_options['MaxCornerAngle']=options.getfieldvalue('MaxCornerAngle',10.)
 	bamg_options['maxnbv']=options.getfieldvalue('maxnbv',10**6)
 	bamg_options['maxsubdiv']=options.getfieldvalue('maxsubdiv',10.)
@@ -313,6 +408,4 @@
 	bamg_options['power']=options.getfieldvalue('power',1.)
 	bamg_options['splitcorners']=options.getfieldvalue('splitcorners',1)
-	bamg_options['geometricalmetric']=options.getfieldvalue('geometricalmetric',0)
-	bamg_options['random']=options.getfieldvalue('rand',True)
 	bamg_options['verbose']=options.getfieldvalue('verbose',1)
 	#}}}
@@ -322,21 +415,62 @@
 
 	# plug results onto model
+	if options.getfieldvalue('vertical',0):
+		md.mesh=mesh2dvertical()
+		md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
+		md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
+		md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+		md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+		md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+		md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
+
+		#Fill in rest of fields:
+		md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
+		md.mesh.numberofvertices=np.size(md.mesh.x)
+		md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
+		md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
+		md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+
+		#Now, build the connectivity tables for this mesh. Doubled in matlab for some reason
+		md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,)
+		md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=1
+
+	elif options.getfieldvalue('3dsurface',0):
+		md.mesh=mesh3dsurface()
+		md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
+		md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
+		md.mesh.z=md.mesh.x
+		md.mesh.z[:]=0
+		md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+		md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+		md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+		md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
+
+		#Fill in rest of fields:
+		md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
+		md.mesh.numberofvertices=np.size(md.mesh.x)
+		md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
+		md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
+		md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+
+	else:
+		md.mesh=mesh2d()
+		md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
+		md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
+		md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+		md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+		md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+		md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
+
+		#Fill in rest of fields:
+		md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
+		md.mesh.numberofvertices=np.size(md.mesh.x)
+		md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
+		md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
+		md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+
+	#Bamg private fields
 	md.private.bamg=OrderedDict()
 	md.private.bamg['mesh']=bamgmesh(bamgmesh_out)
 	md.private.bamg['geometry']=bamggeom(bamggeom_out)
-	md.mesh = mesh2d()
-	md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
-	md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
-	md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
-	md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
-	md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
-	md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
-
-	#Fill in rest of fields:
-	md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
-	md.mesh.numberofvertices=np.size(md.mesh.x)
-	md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
-	md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
-	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
 	md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity
 	md.mesh.elementconnectivity[np.nonzero(np.isnan(md.mesh.elementconnectivity))]=0
@@ -345,5 +479,5 @@
 	#Check for orphan
 	if np.any(np.logical_not(np.in1d(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements.flat))):
-		raise RuntimeError("Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)")
+		raise RuntimeError("Output mesh has orphans. Check your Domain and/or RequiredVertices")
 
 	return md
