demo_qmeshcut_ex1

 qmeshcut demonstration

 by Qianqian Fang, <fangq at nmr.mgh.harvard.edu>

 to demonstrate how to use qmeshcut to produce cross-sectional plot 
 of an un-structured (tetrahedral) mesh

0001 % qmeshcut demonstration
0002 %
0003 % by Qianqian Fang, <fangq at nmr.mgh.harvard.edu>
0004 %
0005 % to demonstrate how to use qmeshcut to produce cross-sectional plot
0006 % of an un-structured (tetrahedral) mesh
0007 
0008 % run vol2mesh demo 1 to create a 3d mesh
0009 
0010 demo_vol2mesh_ex1
0011 
0012 % define a plane by 3 points, in this case, z=mean(node(:,3))
0013 
0014 z0=mean(node(:,3));
0015 
0016 plane=[min(node(:,1)) min(node(:,2)) z0
0017        min(node(:,1)) max(node(:,2)) z0
0018        max(node(:,1)) min(node(:,2)) z0];
0019 
0020 % run qmeshcut to get the cross-section information at z=mean(node(:,1))
0021 % use the x-coordinates as the nodal values
0022 
0023 [cutpos,cutvalue,facedata]=qmeshcut(elem(:,1:4),node,node(:,1),plane);
0024 
0025 % plot your results
0026 
0027 figure;
0028 hsurf=trimesh(face(:,1:3),node(:,1),node(:,2),node(:,3),'facecolor','none');
0029 hold on;
0030 if(isoctavemesh)
0031   hcut=patch('Faces',facedata,'Vertices',cutpos);
0032 else
0033   hcut=patch('Faces',facedata,'Vertices',cutpos,'FaceVertexCData',cutvalue,'facecolor','interp');
0034 end
0035 %set(hcut, 'linestyle','none')
0036 axis equal;
0037 
0038 % qmeshcut can also cut a surface
0039 
0040 [bcutpos,bcutvalue,bcutedges]=qmeshcut(face(:,1:3),node,node(:,1),plane);
0041 [bcutpos,bcutedges]=removedupnodes(bcutpos,bcutedges);
0042 bcutloop=extractloops(bcutedges);
0043 
0044 bcutloop(isnan(bcutloop))=[]; % there can be multiple loops, remove the separators
0045 
0046 % plot the plane-surface cuts
0047 
0048 plot3(bcutpos(bcutloop,1),bcutpos(bcutloop,2),bcutpos(bcutloop,3),'r','LineWidth',4);
0049 
0050 % essencially, this should be the same as you do a removedupnodes(cutpos,facedata)
0051 % and then call extractloop(facedata)
0052 
0053 
0054 % qmeshcut can also cut along an isosurface
0055 
0056 % define a field over the mesh: sensitivity map from a source/detector pair
0057 
0058 r1=[node(:,1)-20,node(:,2)-25,node(:,3)-25];
0059 r2=[node(:,1)-10,node(:,2)-25,node(:,3)-14];
0060 r1=sqrt(r1(:,1).^2+r1(:,2).^2+r1(:,3).^2);
0061 r2=sqrt(r2(:,1).^2+r2(:,2).^2+r2(:,3).^2);
0062 
0063 k=10;
0064 g1=exp(sqrt(-1)*k*r1)./(4*pi*r1); % calculate the Green's function
0065 g2=exp(sqrt(-1)*k*r2)./(4*pi*r2);
0066 g12=g1.*g2;  % this is the sensitivity map
0067 
0068 figure
0069 plotmesh([node log10(abs(g12))],elem,'facealpha',0.5,'linestyle','none'); % plot the mesh
0070 
0071 hold on;
0072 % cut the mesh at value=-4
0073 [cutpos,cutvalue,facedata]=qmeshcut(elem(:,1:4),node(:,1:3),log10(abs(g12)),-4); 
0074 patch('Vertices',cutpos,'Faces',facedata,'FaceVertexCData',cutvalue,'FaceColor','interp');
0075 
0076 % cut the mesh at value=-4.5
0077 [cutpos,cutvalue,facedata]=qmeshcut(elem(:,1:4),node(:,1:3),log10(abs(g12)),-4.5);
0078 patch('Vertices',cutpos,'Faces',facedata,'FaceVertexCData',cutvalue,'FaceColor','interp');
0079