Index: /issm/trunk/src/m/classes/public/plot/plot_segmentonneumann.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_segmentonneumann.m (revision 1615) +++ /issm/trunk/src/m/classes/public/plot/plot_segmentonneumann.m (revision 1616) @@ -1,3 +1,3 @@ -function plot_segmentonneumann(md,options_structure,width,i,datai); +function plot_segmentonneumann(md,options_structure,width,i,data); %PLOT_SEGMENTONNEUMANN - plot segment on neumann BC % @@ -10,89 +10,68 @@ subplot(width,width,i); -if ~isnan(options_structure.unitmultiplier), - md.x=md.x*options_structure.unitmultiplier; - md.y=md.y*options_structure.unitmultiplier; - md.z=md.z*options_structure.unitmultiplier; -end - -if strcmp(datai,'segmentonneumann_diag') segmentonneumann=md.segmentonneumann_diag;end -if strcmp(datai,'segmentonneumann_prog') segmentonneumann=md.segmentonneumann_prog;end - -%we are dealing with the loads. -length_icefront=sqrt( (md.x(segmentonneumann(:,1))-md.x(segmentonneumann(:,2))).^2 + (md.y(segmentonneumann(:,1))-md.y(segmentonneumann(:,2))).^2 ); -normal_icefront(:,1)=cos(atan2( (md.x(segmentonneumann(:,1))-md.x(segmentonneumann(:,2))) , (md.y(segmentonneumann(:,2))-md.y(segmentonneumann(:,1))) ) ); -normal_icefront(:,2)=sin(atan2( (md.x(segmentonneumann(:,1))-md.x(segmentonneumann(:,2))) , (md.y(segmentonneumann(:,2))-md.y(segmentonneumann(:,1))) ) ); +%process mesh and data +[x y z elements is2d]=processmesh(md,options_structure); +if strcmp(data,'segmentonneumann_diag') segmentonneumann=md.segmentonneumann_diag;end +if strcmp(data,'segmentonneumann_prog') segmentonneumann=md.segmentonneumann_prog;end if strcmpi(md.type,'2d'), - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black'); + %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; - pos=segmentonneumann(:,3); - A=md.elements(pos,1); B=md.elements(pos,2); C=md.elements(pos,3); - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black'); + %highlight elements on neumann in Green + pos=segmentonneumann(:,end); + 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','green','EdgeColor','black'); hold on; - x=md.x(md.segmentonneumann_diag(:,1:2))*[1;1]/2; - y=md.y(md.segmentonneumann_diag(:,1:2))*[1;1]/2; + %display arrows pointing outward + xstart=mean(x(segmentonneumann(:,1:end-1)),2); + ystart=mean(y(segmentonneumann(:,1:end-1)),2); + length=sqrt((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))).^2 + (y(segmentonneumann(:,1))-y(segmentonneumann(:,2))).^2 ); + normal(:,1)=cos(atan2((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))) , (y(segmentonneumann(:,2))-y(segmentonneumann(:,1))))); + normal(:,2)=sin(atan2((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))) , (y(segmentonneumann(:,2))-y(segmentonneumann(:,1))))); + xend=xstart+length.*normal(:,1); + yend=ystart+length.*normal(:,2); + q=quiver(xstart,ystart,xend-xstart,yend-ystart); hold on; + h3=plot(xstart,ystart,'r*'); +else - length=sqrt( (md.x(md.segmentonneumann_diag(:,1))-md.x(md.segmentonneumann_diag(:,2))).^2 + (md.y(md.segmentonneumann_diag(:,1))-md.y(md.segmentonneumann_diag(:,2))).^2 ); - normal(:,1)=cos(atan2( (md.x(md.segmentonneumann_diag(:,1))-md.x(md.segmentonneumann_diag(:,2))) , (md.y(md.segmentonneumann_diag(:,2))-md.y(md.segmentonneumann_diag(:,1))) ) ); - normal(:,2)=sin(atan2( (md.x(md.segmentonneumann_diag(:,1))-md.x(md.segmentonneumann_diag(:,2))) , (md.y(md.segmentonneumann_diag(:,2))-md.y(md.segmentonneumann_diag(:,1))) ) ); - - xend=x+length.*normal(:,1); - yend=y+length.*normal(:,2); - quiver(x,y,xend-x,yend-y); hold on; - plot(x,y,'r*'); - legend('element edges on ice front','grids on ice front','normal vectors') -else - for n=1:size(segmentonneumann,1), - hold on + %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; - %Plot the two triangles of each quad (element edges on ice front) - A=segmentonneumann(n,1); B=segmentonneumann(n,2); C=segmentonneumann(n,3); D=segmentonneumann(n,4); E=md.numberofgrids+1; - coordinates= [md.x md.y md.z; mean(md.x(segmentonneumann(n,1:4))) mean(md.y(segmentonneumann(n,1:4))) mean(md.z(segmentonneumann(n,1:4)))]; - p1=patch( 'Faces', [A B E A], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2); - patch( 'Faces', [A D E A], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2); - patch( 'Faces', [D C E D], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2); - patch( 'Faces', [C B E C], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2); + %highlight elements on neumann in Green + pos=segmentonneumann(:,end); + 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','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'); + hold on; - %Plot the nodes of each quad (grids on ice front) - p2=plot3(md.x(segmentonneumann(n,1:4)),md.y(segmentonneumann(n,1:4)),md.z(segmentonneumann(n,1:4)),'r.','markersize',12); - plot3(coordinates(end,1),coordinates(end,2),coordinates(end,3),'g.','markersize',12); + %display arrows pointing outward + xstart=mean(x(segmentonneumann(:,1:end-1)),2); + ystart=mean(y(segmentonneumann(:,1:end-1)),2); + zstart=mean(z(segmentonneumann(:,1:end-1)),2); + length=sqrt((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))).^2 + (y(segmentonneumann(:,1))-y(segmentonneumann(:,2))).^2 ); + normal(:,1)=cos(atan2((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))) , (y(segmentonneumann(:,2))-y(segmentonneumann(:,1))))); + normal(:,2)=sin(atan2((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))) , (y(segmentonneumann(:,2))-y(segmentonneumann(:,1))))); + xend=xstart+length.*normal(:,1); + yend=ystart+length.*normal(:,2); + q=quiver3(xstart,ystart,zstart,xend-xstart,yend-ystart,0); hold on; + h3=plot3(xstart,ystart,zstart,'r*'); +end - %Retrieve coordinates of the four nodes of the Quad - xyz1=[md.x(segmentonneumann(n,1));md.y(segmentonneumann(n,1));md.z(segmentonneumann(n,1));]; - xyz2=[md.x(segmentonneumann(n,2));md.y(segmentonneumann(n,2));md.z(segmentonneumann(n,2));]; - xyz3=[md.x(segmentonneumann(n,3));md.y(segmentonneumann(n,3));md.z(segmentonneumann(n,3));]; - xyz4=[md.x(segmentonneumann(n,4));md.y(segmentonneumann(n,4));md.z(segmentonneumann(n,4));]; - xyz5=coordinates(end,:)'; - - %Build the norms of each triangle - V1=cross(xyz5-xyz1,xyz5-xyz2); normal1=1/norm(V1)*V1'; - V2=cross(xyz5-xyz2,xyz5-xyz3); normal2=1/norm(V2)*V2'; - V3=cross(xyz5-xyz3,xyz5-xyz4); normal3=1/norm(V3)*V3'; - V4=cross(xyz5-xyz4,xyz5-xyz1); normal4=1/norm(V4)*V4'; - - xstart1=mean([xyz1(1) xyz2(1) xyz5(1)]); ystart1=mean([xyz1(2) xyz2(2) xyz5(2)]); zstart1=mean([xyz1(3) xyz2(3) xyz5(3)]); - xstart2=mean([xyz2(1) xyz3(1) xyz5(1)]); ystart2=mean([xyz2(2) xyz3(2) xyz5(2)]); zstart2=mean([xyz2(3) xyz3(3) xyz5(3)]); - xstart3=mean([xyz3(1) xyz4(1) xyz5(1)]); ystart3=mean([xyz3(2) xyz4(2) xyz5(2)]); zstart3=mean([xyz3(3) xyz4(3) xyz5(3)]); - xstart4=mean([xyz4(1) xyz1(1) xyz5(1)]); ystart4=mean([xyz4(2) xyz1(2) xyz5(2)]); zstart4=mean([xyz4(3) xyz1(3) xyz5(3)]); - - xend1=xstart1+length_icefront(n)*normal1(1); yend1=ystart1+length_icefront(n)*normal1(2); zend1=zstart1+length_icefront(n)*normal1(3); - xend2=xstart2+length_icefront(n)*normal2(1); yend2=ystart2+length_icefront(n)*normal2(2); zend2=zstart2+length_icefront(n)*normal2(3); - xend3=xstart3+length_icefront(n)*normal3(1); yend3=ystart3+length_icefront(n)*normal3(2); zend3=zstart3+length_icefront(n)*normal3(3); - xend4=xstart4+length_icefront(n)*normal4(1); yend4=ystart4+length_icefront(n)*normal4(2); zend4=zstart4+length_icefront(n)*normal4(3); - - %plot the normals of each triangle (normal vectors) - p3=plot3([xstart1 xend1],[ystart1 yend1],[zstart1 zend1],'r-'); - plot3([xstart2 xend2],[ystart2 yend2],[zstart2 zend2],'r-'); - plot3([xstart3 xend3],[ystart3 yend3],[zstart3 zend3],'r-'); - plot3([xstart4 xend4],[ystart4 yend4],[zstart4 zend4],'r-'); - - end - legend([p1,p2,p3],'element edges on ice front','grids on ice front','normal vectors'); -end +%legend +legend([h2,q],'element on ice front','normal vectors') %apply options