Changeset 25619

Timestamp:

10/01/20 12:57:00 (4 years ago)

Author:

dlcheng

Message:

CHG: bamg.m discontinous track/domain intersection support.

File:

• issm/trunk-jpl/src/m/mesh/bamg.m (modified) (1 diff)

issm/trunk-jpl/src/m/mesh/bamg.m

@@ -364 +364 @@
         %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