Index: /issm/trunk-jpl/src/m/mesh/bamg.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/bamg.m	(revision 22211)
+++ /issm/trunk-jpl/src/m/mesh/bamg.m	(revision 22212)
@@ -5,4 +5,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)
@@ -59,4 +63,5 @@
 
 subdomain_ref = 1;
+hole_ref = 1;
 % Bamg Geometry parameters {{{
 if exist(options,'domain'),
@@ -82,4 +87,47 @@
 	end
 
+	holes = [];
+	if exist(options,'holes'),
+		holesfile=getfieldvalue(options,'holes');
+		if ischar(holesfile),
+			if ~exist(holesfile,'file') error(['bamg error message: file ' holesfile ' not found']); end
+
+			%read holes according to its extension: 
+			[path,name,ext]=fileparts(holesfile);
+			if strcmp(ext,'.exp'),
+				holes=expread(holesfile);
+			elseif strcmp(ext,'.shp'),
+				holes=shpread(holesfile);
+			else
+				error(['bamg error message: file ' holesfile ' format not supported (.shp or .exp)']);
+			end
+		elseif isstruct(holesfile),
+			holes = holesfile;
+		else
+			error('''holes'' type not supported yet');
+		end
+	end
+	subdomains = [];
+	if exist(options,'subdomains'),
+		subdomainsfile=getfieldvalue(options,'subdomains');
+		if ischar(subdomainsfile),
+			if ~exist(subdomainsfile,'file') error(['bamg error message: file ' subdomainsfile ' not found']); end
+
+			%read subdomains according to its extension: 
+			[path,name,ext]=fileparts(subdomainsfile);
+			if strcmp(ext,'.exp'),
+				subdomains=expread(subdomainsfile);
+			elseif strcmp(ext,'.shp'),
+				subdomains=shpread(subdomainsfile);
+			else
+				error(['bamg error message: file ' subdomainsfile ' format not supported (.shp or .exp)']);
+			end
+		elseif isstruct(subdomainsfile),
+			subdomains = subdomainsfile;
+		else
+			error('''subdomains'' type not supported yet');
+		end
+	end
+
 	%Build geometry 
 	count=0;
@@ -99,9 +147,74 @@
 		end
 
+		%Check orientation
+		nods=domain(i).nods-1; %the domain are closed 1=end;
+		test = sum([(domain(i).x(2:nods+1) - domain(i).x(1:nods)).*(domain(i).y(2:nods+1) + domain(i).y(1:nods))]);
+		if (i==1 && test>0) || (i>1 && test<0),
+			disp('At least one contour was not correctly oriented and has been re-oriented');
+			domain(i).x = flipud(domain(i).x); domain(i).y = flipud(domain(i).y);
+		end
+
 		%Add all points to bamg_geometry
-		nods=domain(i).nods-1; %the domain are closed 1=end;
 		bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]];
+		bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1*ones(nods,1)]];
+		if i>1,
+			bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 -subdomain_ref]; subdomain_ref = subdomain_ref+1;
+		else
+			bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 0];
+		end
+
+		%update counter
+		count=count+nods;
+	end
+	for i=1:length(holes),
+
+		%Check that the subdomains is closed
+		if (holes(i).x(1)~=holes(i).x(end) | holes(i).y(1)~=holes(i).y(end)),
+			error('bamg error message: all contours provided in ''holes'' should be closed');
+		end
+
+		%Checks that all holes are INSIDE the principle domain outline
+		flags=ContourToNodes(holes(i).x,holes(i).y,domain(1),0);
+		if any(~flags), error('bamg error message: All holes should be strictly inside the principal domain'); end
+
+		%Check that hole is correctly oriented
+		nods=holes(i).nods-1; %the holes are closed 1=end;
+		if(sum([(holes(i).x(2:nods+1) - holes(i).x(1:nods)).*(holes(i).y(2:nods+1) + holes(i).y(1:nods))]))<0
+			disp('At least one hole was not correctly oriented and has been re-oriented');
+			holes(i).x = flipud(holes(i).x); holes(i).y = flipud(holes(i).y);
+		end
+
+		%Add all points to bamg_geometry
+		bamg_geometry.Vertices=[bamg_geometry.Vertices; [holes(i).x(1:nods) holes(i).y(1:nods) ones(nods,1)]];
 		bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1.*ones(nods,1)]];
-		if i>1, bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 subdomain_ref]; subdomain_ref = subdomain_ref+1; end
+		bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 -hole_ref]; hole_ref = hole_ref+1;
+
+		%update counter
+		count=count+nods;
+	end
+	for i=1:length(subdomains),
+
+		%Check that the subdomains is closed
+		if (subdomains(i).x(1)~=subdomains(i).x(end) | subdomains(i).y(1)~=subdomains(i).y(end)),
+			error('bamg error message: all contours provided in ''subdomains'' should be closed');
+		end
+
+		%Checks that all holes are INSIDE the principle domain outline
+		flags=ContourToNodes(subdomains(i).x,subdomains(i).y,domain(1),0);
+		if any(~flags),
+			error('bamg error message: All holes should be strictly inside the principal domain');
+		end
+
+		%Check that hole is correctly oriented
+		nods=subdomains(i).nods-1; %the subdomains are closed 1=end;
+		if(sum([(subdomains(i).x(2:nods+1) - subdomains(i).x(1:nods)).*(subdomains(i).y(2:nods+1) + subdomains(i).y(1:nods))]))>0
+			disp('At least one subdomain was not correctly oriented and has been re-oriented');
+			subdomains(i).x = flipud(subdomains(i).x); subdomains(i).y = flipud(subdomains(i).y);
+		end
+
+		%Add all points to bamg_geometry
+		bamg_geometry.Vertices=[bamg_geometry.Vertices; [subdomains(i).x(1:nods) subdomains(i).y(1:nods) ones(nods,1)]];
+		bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1.*ones(nods,1)]];
+		bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 subdomain_ref]; subdomain_ref = subdomain_ref+1;
 
 		%update counter