Index: /issm/trunk-jpl/src/m/mesh/bamg.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/bamg.py	(revision 22298)
+++ /issm/trunk-jpl/src/m/mesh/bamg.py	(revision 22299)
@@ -22,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)
@@ -74,10 +78,9 @@
 	bamg_mesh=bamgmesh()
 
-        subdomain_ref = 1
-        hole_ref = 1
+	subdomain_ref = 1
+	hole_ref = 1
 
 	# Bamg Geometry parameters {{{
 	if options.exist('domain'):
-
 		#Check that file exists
 		domainfile=options.getfieldvalue('domain')
@@ -92,10 +95,9 @@
 		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]
@@ -103,22 +105,101 @@
 					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
 			nods=domaini['nods']-1    #the domain are closed 0=end
 			bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods],domaini['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))
-                        if i:
-                            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]))
-
-			# 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]))
-
+			if i:
+				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):