source: issm/oecreview/Archive/15392-16133/ISSM-15611-15612.diff@ 16134

../trunk-jpl/src/m/plot/plot_icefront.m

@@ -11 +11 @@
 
 %process mesh and data
 [x y z elements is2d isplanet]=processmesh(md,[],options);
-icefront=md.diagnostic.icefront;
+ice=(md.mask.icelevelset>0);
+noice=(md.mask.icelevelset<=0);
+zeroice=(md.mask.icelevelset==0);
+elementice=sum(ice(md.mesh.elements),2);
+elementnoice=sum(noice(md.mesh.elements),2);
+elementzeroice=sum(zeroice(md.mesh.elements),2);
+icefront=(elementice & elementnoice) & ~(elementice==2 & elementzeroice);
 
 if (md.mesh.dimension==2),
 
@@ -21 +27 @@
         hold on;
 
         %highlight elements on neumann
-        pos=find(icefront(:,end)==1);
-        pos=icefront(pos,end-1);
+        pos=find(icefront);
         A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
         h2=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','blue','EdgeColor','black');
-        pos=find(icefront(:,end)==0);
-        pos=icefront(pos,end-1);
-        A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-        h3=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
+        hold on;
 
-        pos=find(icefront(:,end)==2);
-        pos=icefront(pos,end-1);
-        A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-        h3bis=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','cyan','EdgeColor','black');
+        %Plot zero icelevelset line
 
-        %display arrows pointing outward
-        xstart=mean(x(icefront(:,1:end-2)),2);
-        ystart=mean(y(icefront(:,1:end-2)),2);
-        length=sqrt((x(icefront(:,1))-x(icefront(:,2))).^2 + (y(icefront(:,1))-y(icefront(:,2))).^2 );
-        normal(:,1)=cos(atan2((x(icefront(:,1))-x(icefront(:,2))) , (y(icefront(:,2))-y(icefront(:,1)))));
-        normal(:,2)=sin(atan2((x(icefront(:,1))-x(icefront(:,2))) , (y(icefront(:,2))-y(icefront(:,1)))));
-        xend=xstart+length.*normal(:,1);
-        yend=ystart+length.*normal(:,2);
-        q=quiver(xstart,ystart,xend-xstart,yend-ystart); hold on;
-        h4=plot(xstart,ystart,'r*');
+
 else
 
         %plot mesh
@@ -57 +47 @@
         hold on;
 
         %highlight elements on neumann
-        pos=find(icefront(:,end)==1);
-        pos=icefront(pos,end-1);
+        pos=find(icefront);
         A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
         h2=patch( 'Faces', [A B C],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','blue','EdgeColor','black');
         patch( 'Faces', [D E F],  'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','blue','EdgeColor','black');
         patch( 'Faces', [A B E D],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','blue','EdgeColor','black');
         patch( 'Faces', [B E F C],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','blue','EdgeColor','black');
         patch( 'Faces', [C A D F],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','blue','EdgeColor','black');
-        pos=find(icefront(:,end)==0);
-        pos=icefront(pos,end-1);
-        A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-        h3=patch( 'Faces', [A B C],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
-        patch( 'Faces', [D E F],  'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
-        patch( 'Faces', [A B E D],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
-        patch( 'Faces', [B E F C],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
-        patch( 'Faces', [C A D F],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
-        pos=find(icefront(:,end)==2);
-        pos=icefront(pos,end-1);
-        A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-        h3bis=patch( 'Faces', [A B C],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','cyan','EdgeColor','black');
-        patch( 'Faces', [D E F],  'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','cyan','EdgeColor','black');
-        patch( 'Faces', [A B E D],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','cyan','EdgeColor','black');
-        patch( 'Faces', [B E F C],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','cyan','EdgeColor','black');
-        patch( 'Faces', [C A D F],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','cyan','EdgeColor','black');
-
-        %display arrows pointing outward
-        xstart=mean(x(icefront(:,1:end-2)),2);
-        ystart=mean(y(icefront(:,1:end-2)),2);
-        zstart=mean(z(icefront(:,1:end-2)),2);
-        length=sqrt((x(icefront(:,1))-x(icefront(:,2))).^2 + (y(icefront(:,1))-y(icefront(:,2))).^2 );
-        normal(:,1)=cos(atan2((x(icefront(:,1))-x(icefront(:,2))) , (y(icefront(:,2))-y(icefront(:,1)))));
-        normal(:,2)=sin(atan2((x(icefront(:,1))-x(icefront(:,2))) , (y(icefront(:,2))-y(icefront(:,1)))));
-        xend=xstart+length.*normal(:,1);
-        yend=ystart+length.*normal(:,2);
-        q=quiver3(xstart,ystart,zstart,xend-xstart,yend-ystart,zeros(numel(xstart),1)); hold on;
-        h4=plot3(xstart,ystart,zstart,'r*');
 end
 
 %legend (disable warnings)
 warning off
-legend([h2,h3,h3bis,q],'element on ice front (Water)','element on ice front (Air)','element on ice front (Ice)','normal vectors')
+legend([h2],'element on ice front')
 warning on
 
 %apply options