source: issm/trunk/src/m/plot/plot_icefront.m@ 13975

function plot_icefront(md,options,width,i,data)
%PLOT_ICEFRONT - plot segment on neumann BC
%
%   Usage:
%      plot_icefront(md,options,width,i);
%
%   See also: PLOTMODEL

%plot mesh boundaries
subplot(width,width,i); 

%process mesh and data
[x y z elements is2d isplanet]=processmesh(md,[],options);
icefront=md.diagnostic.icefront;

if (md.mesh.dimension==2),

        %plot mesh
        A=elements(:,1); B=elements(:,2); C=elements(:,3); 
        h1=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
        hold on;

        %highlight elements on neumann
        pos=find(icefront(:,end)==1);
        pos=icefront(pos,end-1);
        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');

        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');

        %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
        A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6);
        h1=patch( 'Faces', [A B C],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
        patch( 'Faces', [D E F],  'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
        patch( 'Faces', [A B E D],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
        patch( 'Faces', [B E F C],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
        patch( 'Faces', [C A D F],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
        hold on;

        %highlight elements on neumann
        pos=find(icefront(:,end)==1);
        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);
        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')
warning on

%apply options
options=addfielddefault(options,'title','Neumann boundary conditions');
options=addfielddefault(options,'colorbar',0);
applyoptions(md,[],options);