Changeset 1126

Timestamp:

06/26/09 11:57:18 (16 years ago)

Author:

Eric.Larour

Message:

Greatly accelerated geography using connectivity tables

File:

• issm/trunk/src/m/classes/public/geography2.m (modified) (8 diffs)

issm/trunk/src/m/classes/public/geography2.m

@@ -1 +1 @@
 function md=geography2(md,landname,iceshelfname,icesheetname)
-
 %Get assigned fields
 x=md.x;
@@ -7 +6 @@
 
 %recover elements and grids on land.
-[gridonland,elementonland]=ContourToMesh(elements,x,y,expread(landname,1),'element and node',1);
+[gridonland,elementonland]=ContourToMesh(elements,x,y,expread(landname,1),'element and node',2);
 
 %use segments to make sure gridonland is correct at the boundary:
-gridonland(md.segments(:,1))=1;
-gridonland(md.segments(:,2))=1;
-
+pos=find(md.segmentmarkers==2);
+gridonland(md.segments(pos,1))=1;
+gridonland(md.segments(pos,2))=1;
 %now, any elements with 3 grids on land should be on land:
 elementsonwater=find(~elementonland);
@@ -42 +41 @@
         elementoniceshelf=ones(md.numberofelements,1);
 else
-        [gridoniceshelf,elementoniceshelf]=ContourToMesh(elements,x,y,expread(iceshelfname,1),'element and node',1);
+        [gridoniceshelf,elementoniceshelf]=ContourToMesh(elements,x,y,expread(iceshelfname,1),'element and node',2);
 end
 
@@ -49 +48 @@
         elementonicesheet=zeros(md.numberofelements,1);
 else
-        [gridonicesheet,elementonicesheet]=ContourToMesh(elements,x,y,expread(icesheetname,1),'element and node',1);
+        [gridonicesheet,elementonicesheet]=ContourToMesh(elements,x,y,expread(icesheetname,1),'element and node',2);
 end
 
@@ -65 +64 @@
 gridonwater=double(~gridonland);
 elementonwater=double(~elementonland);
+
 
 %correct for islands:
@@ -89 +89 @@
 
 %swap segments so they are all oriented the same way
+maxcon=size(md.nodeconnectivity,2);
 for i=1:length(segments),
-        pair=find(  (md.elements(:,1)==segments(i,1)  | md.elements(:,2)==segments(i,1)  | md.elements(:,3)==segments(i,1) ) & ...
-                    (md.elements(:,1)==segments(i,2)  | md.elements(:,2)==segments(i,2)  | md.elements(:,3)==segments(i,2) )   ...
-                         );
+
+        pair=intersect(md.nodeconnectivity(segments(i,1),1:md.nodeconnectivity(segments(i,1),maxcon)),md.nodeconnectivity(segments(i,2),1:md.nodeconnectivity(segments(i,2),maxcon)));
+
     if md.elementonwater(pair(1)),
                 segments(i,3)=pair(2);
@@ -98 +99 @@
                 segments(i,3)=pair(1);
         end
+        
         ord1=find(segments(i,1)==md.elements(segments(i,3),:));
         ord2=find(segments(i,2)==md.elements(segments(i,3),:));
@@ -108 +110 @@
         end
         segments(i,1:2)=fliplr(segments(i,1:2));
+        
 end
 md.segmentonneumann_diag=segments;
 md.counter=2;
 md.segmentmarkers(:)=1;
-
-