Changeset 1616


Ignore:
Timestamp:
08/10/09 09:28:58 (16 years ago)
Author:
Mathieu Morlighem
Message:

extended plot_segmentonneumann in 3d

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/m/classes/public/plot/plot_segmentonneumann.m

    r1217 r1616  
    1 function plot_segmentonneumann(md,options_structure,width,i,datai);
     1function plot_segmentonneumann(md,options_structure,width,i,data);
    22%PLOT_SEGMENTONNEUMANN - plot segment on neumann BC
    33%
    … …  
    1010subplot(width,width,i);
    1111
    12 if ~isnan(options_structure.unitmultiplier),
    13         md.x=md.x*options_structure.unitmultiplier;
    14         md.y=md.y*options_structure.unitmultiplier;
    15         md.z=md.z*options_structure.unitmultiplier;
    16 end
    17 
    18 if strcmp(datai,'segmentonneumann_diag') segmentonneumann=md.segmentonneumann_diag;end
    19 if strcmp(datai,'segmentonneumann_prog') segmentonneumann=md.segmentonneumann_prog;end
    20 
    21 %we are dealing with the loads.
    22 length_icefront=sqrt( (md.x(segmentonneumann(:,1))-md.x(segmentonneumann(:,2))).^2 + (md.y(segmentonneumann(:,1))-md.y(segmentonneumann(:,2))).^2 );
    23 normal_icefront(:,1)=cos(atan2( (md.x(segmentonneumann(:,1))-md.x(segmentonneumann(:,2))) , (md.y(segmentonneumann(:,2))-md.y(segmentonneumann(:,1))) ) );
    24 normal_icefront(:,2)=sin(atan2( (md.x(segmentonneumann(:,1))-md.x(segmentonneumann(:,2))) , (md.y(segmentonneumann(:,2))-md.y(segmentonneumann(:,1))) ) );
     12%process mesh and data
     13[x y z elements is2d]=processmesh(md,options_structure);
     14if strcmp(data,'segmentonneumann_diag') segmentonneumann=md.segmentonneumann_diag;end
     15if strcmp(data,'segmentonneumann_prog') segmentonneumann=md.segmentonneumann_prog;end
    2516
    2617if strcmpi(md.type,'2d'),
    2718
    28         A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3);
    29         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
     19        %plot mesh
     20        A=elements(:,1); B=elements(:,2); C=elements(:,3);
     21        h1=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
    3022        hold on;
    3123
    32         pos=segmentonneumann(:,3);
    33         A=md.elements(pos,1); B=md.elements(pos,2); C=md.elements(pos,3);
    34         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
     24        %highlight elements on neumann in Green
     25        pos=segmentonneumann(:,end);
     26        A=elements(pos,1); B=elements(pos,2); C=elements(pos,3);
     27        h2=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
    3528        hold on;
    3629
    37         x=md.x(md.segmentonneumann_diag(:,1:2))*[1;1]/2;
    38         y=md.y(md.segmentonneumann_diag(:,1:2))*[1;1]/2;
     30        %display arrows pointing outward
     31        xstart=mean(x(segmentonneumann(:,1:end-1)),2);
     32        ystart=mean(y(segmentonneumann(:,1:end-1)),2);
     33        length=sqrt((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))).^2 + (y(segmentonneumann(:,1))-y(segmentonneumann(:,2))).^2 );
     34        normal(:,1)=cos(atan2((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))) , (y(segmentonneumann(:,2))-y(segmentonneumann(:,1)))));
     35        normal(:,2)=sin(atan2((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))) , (y(segmentonneumann(:,2))-y(segmentonneumann(:,1)))));
     36        xend=xstart+length.*normal(:,1);
     37        yend=ystart+length.*normal(:,2);
     38        q=quiver(xstart,ystart,xend-xstart,yend-ystart); hold on;
     39        h3=plot(xstart,ystart,'r*');
     40else
    3941
    40         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 );
    41         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))) ) );
    42         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))) ) );
    43          
    44         xend=x+length.*normal(:,1);
    45         yend=y+length.*normal(:,2);
    46         quiver(x,y,xend-x,yend-y); hold on;
    47         plot(x,y,'r*');
    48         legend('element edges on ice front','grids on ice front','normal vectors')
    49 else
    50         for n=1:size(segmentonneumann,1),
    51                 hold on
     42        %plot mesh
     43        A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6);
     44        h1=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
     45        patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
     46        patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
     47        patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
     48        patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor','black');
     49        hold on;
    5250
    53                 %Plot the two triangles of each quad (element edges on ice front)
    54                 A=segmentonneumann(n,1);  B=segmentonneumann(n,2); C=segmentonneumann(n,3); D=segmentonneumann(n,4); E=md.numberofgrids+1;
    55                 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)))];
    56                 p1=patch( 'Faces', [A B E A], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2);
    57                 patch( 'Faces', [A D E A], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2);
    58                 patch( 'Faces', [D C E D], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2);
    59                 patch( 'Faces', [C B E C], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2);
     51        %highlight elements on neumann in Green
     52        pos=segmentonneumann(:,end);
     53        A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
     54        h2=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
     55        patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
     56        patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
     57        patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
     58        patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','green','EdgeColor','black');
     59        hold on;
    6060
    61                 %Plot the nodes of each quad (grids on ice front)
    62                 p2=plot3(md.x(segmentonneumann(n,1:4)),md.y(segmentonneumann(n,1:4)),md.z(segmentonneumann(n,1:4)),'r.','markersize',12);
    63                 plot3(coordinates(end,1),coordinates(end,2),coordinates(end,3),'g.','markersize',12);
     61        %display arrows pointing outward
     62        xstart=mean(x(segmentonneumann(:,1:end-1)),2);
     63        ystart=mean(y(segmentonneumann(:,1:end-1)),2);
     64        zstart=mean(z(segmentonneumann(:,1:end-1)),2);
     65        length=sqrt((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))).^2 + (y(segmentonneumann(:,1))-y(segmentonneumann(:,2))).^2 );
     66        normal(:,1)=cos(atan2((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))) , (y(segmentonneumann(:,2))-y(segmentonneumann(:,1)))));
     67        normal(:,2)=sin(atan2((x(segmentonneumann(:,1))-x(segmentonneumann(:,2))) , (y(segmentonneumann(:,2))-y(segmentonneumann(:,1)))));
     68        xend=xstart+length.*normal(:,1);
     69        yend=ystart+length.*normal(:,2);
     70        q=quiver3(xstart,ystart,zstart,xend-xstart,yend-ystart,0); hold on;
     71        h3=plot3(xstart,ystart,zstart,'r*');
     72end
    6473
    65                 %Retrieve coordinates of the four nodes of the Quad
    66                 xyz1=[md.x(segmentonneumann(n,1));md.y(segmentonneumann(n,1));md.z(segmentonneumann(n,1));];
    67                 xyz2=[md.x(segmentonneumann(n,2));md.y(segmentonneumann(n,2));md.z(segmentonneumann(n,2));];
    68                 xyz3=[md.x(segmentonneumann(n,3));md.y(segmentonneumann(n,3));md.z(segmentonneumann(n,3));];
    69                 xyz4=[md.x(segmentonneumann(n,4));md.y(segmentonneumann(n,4));md.z(segmentonneumann(n,4));];
    70                 xyz5=coordinates(end,:)';
    71                
    72                 %Build the norms of each triangle
    73                 V1=cross(xyz5-xyz1,xyz5-xyz2); normal1=1/norm(V1)*V1';
    74                 V2=cross(xyz5-xyz2,xyz5-xyz3); normal2=1/norm(V2)*V2';
    75                 V3=cross(xyz5-xyz3,xyz5-xyz4); normal3=1/norm(V3)*V3';
    76                 V4=cross(xyz5-xyz4,xyz5-xyz1); normal4=1/norm(V4)*V4';
    77 
    78                 xstart1=mean([xyz1(1) xyz2(1) xyz5(1)]); ystart1=mean([xyz1(2) xyz2(2) xyz5(2)]); zstart1=mean([xyz1(3) xyz2(3) xyz5(3)]);
    79                 xstart2=mean([xyz2(1) xyz3(1) xyz5(1)]); ystart2=mean([xyz2(2) xyz3(2) xyz5(2)]); zstart2=mean([xyz2(3) xyz3(3) xyz5(3)]);
    80                 xstart3=mean([xyz3(1) xyz4(1) xyz5(1)]); ystart3=mean([xyz3(2) xyz4(2) xyz5(2)]); zstart3=mean([xyz3(3) xyz4(3) xyz5(3)]);
    81                 xstart4=mean([xyz4(1) xyz1(1) xyz5(1)]); ystart4=mean([xyz4(2) xyz1(2) xyz5(2)]); zstart4=mean([xyz4(3) xyz1(3) xyz5(3)]);
    82 
    83                 xend1=xstart1+length_icefront(n)*normal1(1);    yend1=ystart1+length_icefront(n)*normal1(2);    zend1=zstart1+length_icefront(n)*normal1(3);
    84                 xend2=xstart2+length_icefront(n)*normal2(1);    yend2=ystart2+length_icefront(n)*normal2(2);    zend2=zstart2+length_icefront(n)*normal2(3);
    85                 xend3=xstart3+length_icefront(n)*normal3(1);    yend3=ystart3+length_icefront(n)*normal3(2);    zend3=zstart3+length_icefront(n)*normal3(3);
    86                 xend4=xstart4+length_icefront(n)*normal4(1);    yend4=ystart4+length_icefront(n)*normal4(2);    zend4=zstart4+length_icefront(n)*normal4(3);
    87                
    88                 %plot the normals of each triangle (normal vectors)
    89                 p3=plot3([xstart1 xend1],[ystart1 yend1],[zstart1 zend1],'r-');
    90                 plot3([xstart2 xend2],[ystart2 yend2],[zstart2 zend2],'r-');
    91                 plot3([xstart3 xend3],[ystart3 yend3],[zstart3 zend3],'r-');
    92                 plot3([xstart4 xend4],[ystart4 yend4],[zstart4 zend4],'r-');
    93 
    94         end
    95         legend([p1,p2,p3],'element edges on ice front','grids on ice front','normal vectors');
    96 end
     74%legend
     75legend([h2,q],'element on ice front','normal vectors')
    9776
    9877%apply options
Note: See TracChangeset for help on using the changeset viewer.