Index: issm/trunk-jpl/src/m/mesh/bamg.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/bamg.m	(revision 24865)
+++ issm/trunk-jpl/src/m/mesh/bamg.m	(revision 25455)
@@ -80,5 +80,5 @@
 			domain=shpread(domainfile);
 		else
-			error(['bamg error message: file ' domainfile ' format not supported (.shp or .exp)']);
+			error(['bamg error message: file ' domainfile ' format not supported (.exp or .shp)']);
 		end
 	elseif isstruct(domainfile),
@@ -140,5 +140,5 @@
 		end
 
-		%Checks that all holes are INSIDE the principle domain outline
+		%Check that all holes are INSIDE the principle domain outline
 		if i>1,
 			flags=ContourToNodes(domain(i).x,domain(i).y,domain(1),0);
@@ -149,5 +149,5 @@
 
 		%Check orientation
-		nods=domain(i).nods-1; %the domain are closed 1=end;
+		nods=domain(i).nods-1; %the domain is closed (domain[1] = domain[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),
@@ -163,11 +163,13 @@
 		bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1*ones(nods,1)]];
 
-		%Flag how many edges we have now, that way we know which edges belong to the domain, will 
-		%be used later for required edges when NoBoundaryRefinement is one: 
+		% Flag how many edges we have now, that way we know which edges belong 
+		% to the domain. Will be used later for required edges if 
+		% NoBoundaryRefinement equals 1.
 		new_edge_length=length(bamg_geometry.Edges);
 		edges_required=(edge_length+1):new_edge_length;
 
 		if i>1,
-			bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 -subdomain_ref]; subdomain_ref = subdomain_ref+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];
@@ -179,15 +181,15 @@
 	for i=1:length(holes),
 
-		%Check that the subdomains is closed
+		%Check that the hole 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
+		%Check that all holes are INSIDE the principal domain (principal domain should be index 0)
 		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;
+		nods=holes(i).nods-1; %the hole is closed (hole[1] = hole[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');
@@ -205,20 +207,21 @@
 	for i=1:length(subdomains),
 
-		%Check that the subdomains is closed
+		%Check that the subdomain 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
+		%Checks that all subdomains are INSIDE the principal domain (principal domain should be index 0)
 		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;
+			error('bamg error message: All subdomains should be strictly inside the principal domain');
+		end
+
+		%Check that subdomain is correctly oriented
+		nods=subdomains(i).nods-1; % the subdomains are closed (subdomains[1] = subdomains[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);
+			subdomains(i).x = flipud(subdomains(i).x); 
+			subdomains(i).y = flipud(subdomains(i).y);
 		end
 
@@ -227,9 +230,11 @@
 		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;
+		bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 subdomain_ref];
+		subdomain_ref = subdomain_ref+1;
 		
 		%update counter
 		count=count+nods;
 	end
+
 	if getfieldvalue(options,'vertical',0),
 		if numel(getfieldvalue(options,'Markers',[]))~=size(bamg_geometry.Edges,1),
@@ -243,16 +248,6 @@
 	%take care of rifts
 	if exist(options,'rifts'),
-
-		%Check that file exists
+		%read rift file according to its extension: 
 		riftfile=getfieldvalue(options,'rifts');
-		[pathr,namer,extr]=fileparts(riftfile);
-		if ~exist(riftfile,'file')
-			error(['bamg error message: file ' riftfile ' not found ']);
-		elseif strcmp(extr,'.exp'),
-			rift=expread(riftfile);
-		elseif strcmp(extr,'.shp'),
-			rift=shpread(riftfile);
-		end
-		%read rift file according to its extension: 
 		[path,name,ext]=fileparts(riftfile);
 		if strcmp(ext,'.exp'),
@@ -261,5 +256,5 @@
 			rift=shpread(riftfile);
 		else
-			error(['bamg error message: file ' riftfile ' format not supported (.shp or .exp)']);
+			error(['bamg error message: file ' riftfile ' format not supported (.exp or .shp)']);
 		end
 
@@ -272,10 +267,10 @@
 
 			elseif any(~flags),
-				%We LOTS of work to do
+				%We have LOTS of work to do
 				disp('Rift tip outside of or on the domain has been detected and is being processed...');
 
 				%check that only one point is outside (for now)
 				if sum(~flags)~=1,
-					error('bamg error message: only one point outside of the domain is supported yet');
+					error('bamg error message: only one point outside of the domain is supported at this time');
 				end
 
@@ -290,7 +285,9 @@
 				end
 
-				%Get cordinate of intersection point
-				x1=rift(i).x(1); y1=rift(i).y(1);
-				x2=rift(i).x(2); y2=rift(i).y(2);
+				%Get coordinate of intersection point
+				x1=rift(i).x(1);
+				y1=rift(i).y(1);
+				x2=rift(i).x(2);
+				y2=rift(i).y(2);
 				for j=1:length(domain(1).x)-1;
 					if SegIntersect([x1 y1; x2 y2],[domain(1).x(j) domain(1).y(j); domain(1).x(j+1) domain(1).y(j+1)]),
@@ -489,8 +486,6 @@
 	md.mesh.numberofvertices=length(md.mesh.x);
 	md.mesh.numberofedges=size(md.mesh.edges,1);
-	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
-
-	%Now, build the connectivity tables for this mesh.
-	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
 
 elseif getfieldvalue(options,'3dsurface',0),
@@ -509,5 +504,6 @@
 	md.mesh.numberofvertices=length(md.mesh.x);
 	md.mesh.numberofedges=size(md.mesh.edges,1);
-	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
 
 else 
@@ -524,5 +520,6 @@
 	md.mesh.numberofvertices=length(md.mesh.x);
 	md.mesh.numberofedges=size(md.mesh.edges,1);
-	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
 end
 
