Index: /issm/trunk-jpl/src/m/mesh/bamg.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/bamg.m	(revision 25618)
+++ /issm/trunk-jpl/src/m/mesh/bamg.m	(revision 25619)
@@ -364,26 +364,52 @@
 	%Deal with tracks
 	if exist(options,'tracks'),
-
 		%read tracks
-		track=getfieldvalue(options,'tracks');
-		if all(ischar(track)),
-			A=expread(track);
-			track=[];
-			for i=1:length(A), 
-				track=[track; [A(i).x A(i).y]];
-			end
-		else
-			track=double(track); %for some reason, it is of class "single"
-		end
+		intrack=getfieldvalue(options,'tracks');
+		if all(ischar(intrack)),
+			intrack=expread(intrack);
+        else
+            intrack=double(intrack); %for some reason, it is of class "single"
+            intrack=struct('x',intrack(:,1),'y',intrack(:,2));
+            intrack=repmat(intrack, [1,1]);
+        end
+        
+        %Process multiple tracks and handle domain intersections
+        track=[];
+        edges=[];
+        newcount=count;
+        for i=1:length(intrack), 
+            nods=size(intrack(i).x,1);
+            newtrack=[intrack(i).x intrack(i).y];
+            newedges=[transpose(newcount+1:newcount+nods-1) transpose(newcount+2:newcount+nods) 3.*ones(nods-1,1)];
+
+            %only keep those inside
+            flags=ContourToNodes(newtrack(:,1),newtrack(:,2),domainfile,0);
+            edgeflags = flags(1:end-1) & flags(2:end);
+
+            %calculate edge offset, accounting for broken tracks
+            offsets=zeros(size(edgeflags,1)+1,3);
+            for j=2:size(offsets,1),
+                if edgeflags(j-1)==0,
+                    offsets(j,:)=offsets(j-1,:)+[1 1 0];
+                else
+                    offsets(j,:)=offsets(j-1,:);
+                end
+            end
+            newedges=newedges-offsets(2:end,:);
+
+            %add track segments
+            newtrack=newtrack(find(flags),:);
+            newedges=newedges(find(edgeflags),:);
+            nods=size(newtrack,1);
+            track=[track; newtrack];
+            edges=[edges; newedges];
+            newcount=newcount+nods;
+        end
 		if(size(track,2)==2), track=[track 3.*ones(size(track,1),1)]; end
-
-		%only keep those inside
-		flags=ContourToNodes(track(:,1),track(:,2),domainfile,0);
-		track=track(find(flags),:);
-
+        
 		%Add all points to bamg_geometry
 		nods=size(track,1);
 		bamg_geometry.Vertices=[bamg_geometry.Vertices; track];
-		bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods-1) transpose(count+2:count+nods)  3.*ones(nods-1,1)]];
+		bamg_geometry.Edges=[bamg_geometry.Edges; edges];
 
 		%update counter