Changeset 3116

Timestamp:

02/24/10 11:30:51 (15 years ago)

Author:

seroussi

Message:

improved geography2

File:

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

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

@@ -25 +25 @@
 else
         error('Invalid area option option');
+end
+
+%Now, build the connectivity tables for this mesh.
+if size(md.nodeconnectivity,1)~=md.numberofgrids,
+        md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
+end
+if size(md.elementconnectivity,1)~=md.numberofelements,
+        md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
 end
 
@@ -87 +95 @@
 elementonwater=double(~elementonland);
 
-
 %correct for islands:
 gridoniceshelf=double(gridoniceshelf & ~gridonicesheet);
@@ -96 +103 @@
 elementonicesheet=double(~elementoniceshelf & ~elementonwater);
 
+%Deal with segments on neumann:
 
+%Get current connectivity
+elementconnectivity=md.elementconnectivity;
 
-%Deal with segments on neumann: 
-elementconnectivity=md.elementconnectivity;
+%put 0 for elements on water
 pos=find(elementconnectivity);
-elementconnectivity(pos)=elementonwater(elementconnectivity(pos));
+elementconnectivity(pos)=elementconnectivity(pos).*(~elementonwater(elementconnectivity(pos)));
 
+%put line of ones for elements on water
 pos=find(elementonwater);
-elementconnectivity(pos,:)=0;
+elementconnectivity(pos,:)=1;% line of ones for elements on water so they won't be considered
 
-num_segments=sum(sum(elementconnectivity));
-segments=zeros(num_segments,3);
+%resort lines (zeros must be at the last column for findsegments)
+elementconnectivity=sort(elementconnectivity,2,'descend');
 
-icefrontelementsonland=find(sum(elementconnectivity,2));
-
-count=1;
-for i=1:numel(icefrontelementsonland),
-        el1=icefrontelementsonland(i);
-        els2=md.elementconnectivity(el1,find(elementconnectivity(el1,:)));
-        for j=1:numel(els2),
-                el2=els2(j);
-                %we have a segment between land el1 and water el2:
-                segments(count,:)=[intersect(md.elements(el1,:),md.elements(el2,:)) el1];
-                
-                ord1=find(segments(count,1)==md.elements(el1,:));
-                ord2=find(segments(count,2)==md.elements(el1,:));
-
-                %swap segment grids if necessary
-                if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
-                        temp=segments(count,1);
-                        segments(count,1)=segments(count,2);
-                        segments(count,2)=temp;
-                end
-                segments(count,1:2)=fliplr(segments(count,1:2));
-                
-                count=count+1;
-        end
-end
-
+%call findsegments to build segment using THIS conectivity
+md.segments=findsegments(md,'elementconnectivity',elementconnectivity);
 
 %some final checks: