1+ function [hs ,hq ] = cartesianUnstructuredMeshPlotSolution(GeomTRI ,mesh ,IntPlane ,Mx_arr ,My_arr ,Mz_arr ,hF )
2+ %
3+ % CartesianUnstructuredMeshPlotSolution(GeomTRI,mesh,IntPlane,Mx_arr,My_arr,Mz_arr,hF)
4+ %
5+ % Intersects the geometry with a plane, and shows a cutaway view of the
6+ % solution.
7+ %
8+ % GeomTri is the triangulation of the geometry (can be obtained using
9+ % convhulln. For this function to work, the geometry is assumed to be
10+ % convex anyway.
11+ % mesh is the Cartesian Unstructured Mesh corresponding to the solution
12+ % IntPlane is a structure with two fields (each of which is a 1 by 3
13+ % array): IntPlane.Normal is the normal to the plane
14+ % IntPlane.Point is a point passing by the plane
15+ % Mx_arr, My_arr, and Mz_arr are the arrays of the solution
16+ %
17+ % Requires the PDE toolbox (used for creating a tetrahedral mesh -
18+ % and the correponding triangular surface mesh on which
19+ % the solution is plotted)
20+
21+ %%
22+ X = [min(GeomTRI .Points(: ,1 )),max(GeomTRI .Points(: ,1 ))] ;
23+ Y = [min(GeomTRI .Points(: ,2 )),max(GeomTRI .Points(: ,2 ))] ;
24+ Z = [min(GeomTRI .Points(: ,3 )),max(GeomTRI .Points(: ,3 ))] ;
25+ IntPlaneN = norm(IntPlane .Normal ) ;
26+ IntPlane.Normal = IntPlane .Normal ./ IntPlaneN ;
27+
28+ %%
29+ ik = [1 ,2 ,3 ,1 ] ;
30+
31+ % keep only the points on the negative side of the plane
32+ iIn = sum((GeomTRI .Points - repmat(IntPlane .Point ,size(GeomTRI .Points ,1 ),1 )).*repmat(IntPlane .Normal ,size(GeomTRI .Points ,1 ),1 ),2 )<0 ;
33+ % Cycle over each segment
34+ jN = 0 ; clear xNew
35+ for j= 1 : size(GeomTRI .ConnectivityList ,1 )
36+ for k= 1 : 3
37+ % find the intersection point between the segment and the
38+ % intersecting plane
39+ xA = GeomTRI .Points(GeomTRI .ConnectivityList(j ,ik(k )),: ) ;
40+ xB = GeomTRI .Points(GeomTRI .ConnectivityList(j ,ik(k + 1 )),: ) ;
41+ t = - sum((xA - IntPlane .Point ).*IntPlane .Normal )/sum((xB - xA ).*IntPlane .Normal ) ;
42+ xInt = xA + t .*(xB - xA ) ;
43+ if (t >= 0 ) & (t <= 1 ) ; plot3(xInt(1 ),xInt(2 ),xInt(3 ),' .'
44+ jN = jN + 1 ;
45+ xNew(jN ,: ) = xInt ;
46+
47+
48+ end
49+ end
50+ end
51+
52+ BisectGeomTRI.Points = GeomTRI .Points ;
53+ BisectGeomTRI .Points(find(~iIn ),: ) = [] ;
54+ BisectGeomTRI.Points = [BisectGeomTRI .Points ;xNew ] ;
55+
56+ BisectGeomTRI.ConnectivityList = convhulln(BisectGeomTRI .Points ) ;
57+
58+
59+ BiSectVorStruct.pos = sum(BisectGeomTRI .Points ,1 ) ;
60+ BiSectVorStruct.vorvx{1 } = BisectGeomTRI .Points ;
61+ BiSectVorStruct.ThatK{1 } = BisectGeomTRI .ConnectivityList ;
62+
63+ VoronoiStruct.pos = sum(GeomTRI .Points ,1 ) ;
64+ VoronoiStruct.vorvx{1 } = GeomTRI .Points ;
65+ VoronoiStruct.ThatK{1 } = GeomTRI .ConnectivityList ;
66+
67+
68+ % PlotVoronoiGeometry(BiSectVorStruct,X,Y,Z,'BisectedHexagonalPlatelet01',hF,[.5,.5,.5])
69+ BisectGeomTRI = triangulation(BisectGeomTRI .ConnectivityList ,BisectGeomTRI .Points ) ;
70+ stlwrite(BisectGeomTRI ,' BisectGeomTRI.stl'
71+ model = createpde(1 );
72+ importGeometry(model ,' BisectGeomTRI.stl'
73+ delete(' BisectGeomTRI.stl'
74+ generateMesh(model ,' GeometricOrder' ' linear' ' Hmax' X )/50 ) ;
75+ hF0 = figure ;
76+ h = pdeplot3D(model ) ;
77+ SurfMesh = triangulation(h .Faces ,h .Vertices ) ;
78+ delete(hF0 )
79+
80+ %% Mesh for quiver plots
81+ generateMesh(model ,' GeometricOrder' ' linear' ' Hmin' X )/10 ) ;
82+ hF0 = figure ;
83+ h = pdeplot3D(model ) ;
84+ QuiverMesh = triangulation(h .Faces ,h .Vertices ) ;
85+ inQuiv = sum((QuiverMesh .Points - repmat(IntPlane .Point ,size(QuiverMesh .Points ,1 ),1 )).*repmat(IntPlane .Normal ,size(QuiverMesh .Points ,1 ),1 ),2 )> (-.01 * diff(X )) ;
86+ delete(hF0 )
87+ %%
88+
89+ Fmx = scatteredInterpolant(mesh .pos_out(: ,1 ),mesh .pos_out(: ,2 ),mesh .pos_out(: ,3 ),Mx_arr ) ;
90+ Fmy = scatteredInterpolant(mesh .pos_out(: ,1 ),mesh .pos_out(: ,2 ),mesh .pos_out(: ,3 ),My_arr ) ;
91+ Fmz = scatteredInterpolant(mesh .pos_out(: ,1 ),mesh .pos_out(: ,2 ),mesh .pos_out(: ,3 ),Mz_arr ) ;
92+
93+ Mx = Fmx(SurfMesh .Points(: ,1 ),SurfMesh .Points(: ,2 ),SurfMesh .Points(: ,3 )) ;
94+ My = Fmy(SurfMesh .Points(: ,1 ),SurfMesh .Points(: ,2 ),SurfMesh .Points(: ,3 )) ;
95+ Mz = Fmz(SurfMesh .Points(: ,1 ),SurfMesh .Points(: ,2 ),SurfMesh .Points(: ,3 )) ;
96+
97+
98+ MxQ = Fmx(QuiverMesh .Points(: ,1 ),QuiverMesh .Points(: ,2 ),QuiverMesh .Points(: ,3 )) ;
99+ MyQ = Fmy(QuiverMesh .Points(: ,1 ),QuiverMesh .Points(: ,2 ),QuiverMesh .Points(: ,3 )) ;
100+ MzQ = Fmz(QuiverMesh .Points(: ,1 ),QuiverMesh .Points(: ,2 ),QuiverMesh .Points(: ,3 )) ;
101+ MxQ(~inQuiv ) = 0 ; MyQ(~inQuiv ) = 0 ; MzQ(~inQuiv ) = 0 ;
102+ QuiverMesh2.Points = QuiverMesh .Points + (.01 * diff(X )).*repmat(IntPlane .Normal ,size(QuiverMesh .Points ,1 ),1 ) ;
103+ ThisColTriangle = ColorFromHorPsiTheta01(Mx ,My ,Mz ) ;
104+ %%
105+
106+
107+ if ~exist(' hF' ' var'
108+ hF = figure(' position' 0 0 600 600 ],' Color' 1 1 1 ]);
109+ ppsz = 1 .*[20 ,19 ] ;
110+ ppps = 1 .*[0 ,0 ,20 ,19 ] ;
111+ set(gcf ,' PaperUnits' ' centimeters' ' PaperSize' ppsz ,' PaperPosition' ppps ) ;
112+ else
113+ ppsz= get(gcf ,' PaperSize'
114+ ppps= get(gcf ,' PaperPosition'
115+ if isequal(get(hF ,' type' ' figure'
116+ set(hF ,' PaperUnits' ' centimeters'
117+
118+ figure(hF )
119+ % clf(hF)
120+ else
121+ axes(hF ) ;
122+ end
123+ end
124+
125+ hold on
126+
127+ axis(' equal'
128+
129+ %%
130+
131+
132+ % PlotVoronoiGeometry(VoronoiStruct,X,Y,Z,'HexagonalPlatelet01',hF,[.5,.5,.5])
133+ PlotOnlyRealEdges(VoronoiStruct .pos ,VoronoiStruct .vorvx ,VoronoiStruct .ThatK ,1 )
134+ PlotOnlyRealEdges(BiSectVorStruct .pos ,BiSectVorStruct .vorvx ,BiSectVorStruct .ThatK ,2 )
135+
136+
137+ hs = trisurf(SurfMesh .ConnectivityList ,SurfMesh .Points(: ,1 ),SurfMesh .Points(: ,2 ),SurfMesh .Points(: ,3 ),' FaceVertexCData' ThisColTriangle ,' linestyle' ' none' ' facelighting' ' flat' ' SpecularStrength' 0 ,' facealpha' 1 ) ;
138+ hq = quiver3(QuiverMesh2 .Points(: ,1 ),QuiverMesh2 .Points(: ,2 ),QuiverMesh2 .Points(: ,3 ),MxQ ,MyQ ,MzQ ,' color' ' k' ' LineWidth' 1.5 ) ;
139+
140+ view(90 ,30 ) ;
141+ light(' position' 0 ,0 ,2 ].*[X(end ),Y(end ),Z(end )])
142+ light(' position' 0 ,-2 ,0 ].*[X(end ),Y(end ),Z(end )],' color' .3 .*[1 ,1 ,1 ])
143+
144+ set(gca ,' xlim' 1 ),X(end )],...
145+ ' ylim' 1 ),Y(end )],...
146+ ' zlim' 1 ),Z(end )],' visible' ' off'
147+
148+ % xlabel('x') ; ylabel('y') ; zlabel('z') ;
149+
150+ end
151+
152+ function PlotOnlyRealEdges(pos ,vorvx ,ThatK ,LW )
153+
154+ TotFaces = 0 ;
155+ if ~exist(' col' ' var'
156+ col = [hsv(size(pos ,1 ));[1 ,1 ,1 ].*0.4 ];
157+ end
158+ for i = 1 : size(pos ,1 )
159+ K = ThatK{i } ;
160+ % K = convhulln(vorvx{i});
161+ % [K,vorvx{i}] = CleanUp(K,vorvx{i}) ;
162+ TheSharedEdges{i } = zeros(size(K ,1 ),size(K ,1 )) ;
163+ % trisurf(K,vorvx{i}(:,1),vorvx{i}(:,2),vorvx{i}(:,3),'FaceColor',col(i,:),'FaceAlpha',.5,'EdgeAlpha',0,'facelighting','flat','SpecularStrength',0)
164+ hold on ;
165+ % stlwrite(TR,['testSTLwrite',num2str(i),'.stl']) ;
166+ for k= 1 : size(K ,1 ) % cycle over triangles
167+ nn{i }{k } = cross((vorvx{i }(K(k ,2 ),: )-vorvx{i }(K(k ,1 ),: )),(vorvx{i }(K(k ,3 ),: )-vorvx{i }(K(k ,1 ),: ))) ;
168+ nn{i }{k } = nn{i }{k }./norm( nn{i }{k }) ;
169+ % line(sigma.*(vorvx{i}(K(k,:),1)-pos(i,1))+pos(i,1),sigma.*(vorvx{i}(K(k,:),2)-pos(i,2)) + pos(i,2),sigma.*(vorvx{i}(K(k,:),3)-pos(i,3)) + pos(i,3))
170+ end
171+ ThatDot = zeros(size(K ,1 ),size(K ,1 )) ;
172+ for k1= 1 : size(K ,1 )
173+ for k2= (k1 + 1 ): size(K ,1 )
174+ ThatDot(k1 ,k2 ) = dot( nn{i }{k1 }, nn{i }{k2 } ) ;
175+ end
176+ end
177+ ThatDot = ThatDot + ThatDot .' ; % triangles belonging to the same planes
178+ ijk = [1 ,2 ,3 ,1 ] ;
179+ CommonEdges = zeros(size(K ,1 ),3 ) ;
180+ for k1= 1 : size(K ,1 )
181+ for k2= 1 : size(K ,1 )
182+ [C ,ia ,ib ] = intersect( K(k1 ,: ),K(k2 ,: )) ;
183+ if numel(C )==2 & abs(abs(ThatDot(k1 ,k2 ))-1 )<1e- 9
184+ ia = sort(ia ) ;
185+ ib = sort(ib ) ;
186+ TheSide1 = isequal(ia ,[1 ;2]).*1 + isequal(ia ,[2 ;3]).*2 + isequal(ia ,[1 ;3]).*3 ;
187+ TheSide2 = isequal(ib ,[1 ;2]).*1 + isequal(ib ,[2 ;3]).*2 + isequal(ib ,[1 ;3]).*3 ;
188+ CommonEdges(k1 ,TheSide1 ) = 1 ;
189+ CommonEdges(k2 ,TheSide2 ) = 1 ;
190+ TheSharedEdges{i }(k1 ,k2 ) = 1 ; TheSharedEdges{i }(k2 ,k1 ) = 1 ;
191+ % TheSharedEdges{i}(k2,TheSide2) = k1 ;
192+ ' '
193+ end
194+ end
195+ end
196+
197+ ' '
198+ % FindPolygons01(TheSharedEdges{i})
199+ CCC = zeros(max(K(: )),max(K(: ))) ;
200+ for k= 1 : size(K ,1 )
201+ for j= 1 : 3
202+ if ~CommonEdges(k ,j )
203+ CCC(K(k ,ijk(j )),K(k ,ijk(j + 1 ))) = 1 ;
204+ CCC(K(k ,ijk(j + 1 )),K(k ,ijk(j ))) = 1 ;
205+ % line(vorvx{i}(K(k,ijk(j:j+1)),1),vorvx{i}(K(k,ijk(j:j+1)),2),vorvx{i}(K(k,ijk(j:j+1)),3),'color','k','linewidth',2)
206+ end
207+ end
208+ end
209+ bins = conncomp(graph(TheSharedEdges{i })) ;
210+ for kb = 1 : max(bins )
211+ BB = bins == kb ; kB = K(BB ,: ) ;
212+ [ukB ,ikB ,ckB ] = unique(kB(: )) ;
213+
214+ GG = graph(CCC(ukB ,ukB )) ;
215+ nOrd = dfsearch(GG ,1 ) ;
216+ iiord = kB(ikB(nOrd )) ;
217+ iiord = [iiord(: );iiord(1 )];
218+ AllFaces{i ,kb } = iiord ;
219+ TotFaces = TotFaces + 1 ;
220+ AllSingleFaces{TotFaces } = iiord ;
221+ AllSingleFacesZones(TotFaces ) = i ;
222+ line(vorvx{i }(iiord ,1 ),vorvx{i }(iiord ,2 ),vorvx{i }(iiord ,3 ),' color' ' k' ' linewidth' LW )
223+ end
224+ ' '
225+
226+
227+ end
228+ end
0 commit comments