Changeset 1223

Timestamp:

07/02/09 15:50:07 (16 years ago)

Author:

Eric.Larour

Message:

Better detection of icefronts

File:

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

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

@@ -1 +1 @@
 function md=geography2(md,landname,iceshelfname,icesheetname)
+
 %Get assigned fields
 x=md.x;
@@ -8 +9 @@
 [gridonland,elementonland]=ContourToMesh(elements,x,y,expread(landname,1),'element and node',2);
 
-%use segments to make sure gridonland is correct at the boundary:
-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:
+%any element with 3 grids on land should be on land:
 elementsonwater=find(~elementonland);
 wrongelements=elementsonwater(find(( gridonland(md.elements(elementsonwater,1)) + gridonland(md.elements(elementsonwater,2)) + gridonland(md.elements(elementsonwater,3)) ...
@@ -18 +15 @@
 elementonland(wrongelements)=1;
 
-%finally, any element with its barycentre in the water should be in the water!
-xelem=x(md.elements)*[1;1;1]/3; yelem=y(md.elements)*[1;1;1]/3;
+%any element with its barycentre on land should be on land:
+weights={[1;1;1],[2;1;1],[1;2;1],[1;1;2]};
+for i=1:length(weights),
+        xelem=x(md.elements)*weights{i}/sum(weights{i});
+        yelem=y(md.elements)*weights{i}/sum(weights{i});
+end
 baryonland=ContourToNodes(xelem,yelem,expread(landname,1),1);
 pos=find(~baryonland); elementonland(pos)=0;
+pos=find(baryonland); elementonland(pos)=1;
 
 %figure out which elements on land are actually in the middle of the ocean!
@@ -29 +31 @@
 connectedtolandsum=sum(connectedtoland,2);
 waterelements=pos1(find(connectedtolandsum==3));
+elementonland(waterelements)=0;
 
-%now use waterelements to correct elementonland: 
-elementonland(waterelements)=0;
+%figure out which elements on water  are actually in the middle of the land!
+pos1=find(~elementonland); 
+connectedtowater=md.elementconnectivity(pos1,:);
+pos=find(connectedtowater); connectedtowater(pos)=elementonland(connectedtowater(pos));
+connectedtowatersum=sum(connectedtowater,2);
+landelements=pos1(find(connectedtowatersum==3));
+elementonland(landelements)=1;
 
 %recover arrays of ice shelf grids and elements, and ice sheet grids and elements.
@@ -74 +82 @@
 elementonicesheet=double(~elementoniceshelf & ~elementonwater);
 
+
+
+%Deal with segments on neumann: 
+elementconnectivity=md.elementconnectivity;
+pos=find(elementconnectivity);
+elementconnectivity(pos)=elementonwater(elementconnectivity(pos));
+
+pos=find(elementonwater);
+elementconnectivity(pos,:)=0;
+
+num_segments=sum(sum(elementconnectivity));
+segments=zeros(num_segments,3);
+
+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
+
 %Return: 
 md.gridoniceshelf=gridoniceshelf;
@@ -84 +131 @@
 md.elementonicesheet=elementonicesheet;
 
-%Deal with segments on neumann: 
-pos=find(md.segmentmarkers==2); %internal segments.
-segments=md.segments(pos,:);
-
-%swap segments so they are all oriented the same way
-maxcon=size(md.nodeconnectivity,2);
-for i=1:length(segments),
-
-        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);
-        else
-                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),:));
-        
-        %swap segment grids if necessary
-        if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
-                temp=segments(i,1);
-                segments(i,1)=segments(i,2);
-                segments(i,2)=temp;
-        end
-        segments(i,1:2)=fliplr(segments(i,1:2));
-        
-end
 md.segmentonneumann_diag=segments;
 md.counter=2;