source: issm/trunk/src/m/model/mesh/rifts/meshyamsrecreateriftsegments.m@ 9619

Last change on this file since 9619 was 9619, checked in by Mathieu Morlighem, 14 years ago

Added rifts object
File size: 3.1 KB
Line 
1function md=meshyamsrecreateriftsegments(md)
2
3 %recreate rift segments: just used for yams. temporaroy routine.
4 pos_record=[];
5 if md.rifts.numrifts,
6 for i=1:md.rifts.numrifts,
7 rift=md.rifts.riftstruct(i);
8
9 %closed rifts first:
10 if length(rift.tips)==2,
11
12 %find tip1 and tip2 for this rift, in the new mesh created by yams.
13 pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
14 tip1=md.segments(pos,1);
15 pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));
16 tip2=md.segments(pos,1);
17
18 %start from tip1, and build segments of this rift.
19 pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
20 pos_record=[pos_record; pos];
21 riftsegs=md.segments(pos,:);
22 while 1,
23 A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
24 %find other segment that holds B.
25 pos=find(md.segments(:,1)==B);
26 pos_record=[pos_record; pos];
27 riftsegs=[riftsegs; md.segments(pos,:)];
28 if riftsegs(end,2)==tip1,
29 break;
30 end
31 end
32 md.rifts.riftstruct(i).segments=riftsegs;
33 md.rifts.riftstruct(i).tips=[tip1 tip2];
34
35 else
36 %ok, this is a rift that opens up to the domain outline. One tip is going to be
37 %double, the other one, single. We are going to start from the single tip, towards the two
38 %other doubles
39
40 %find tip1 and tip2 for this rift, in the new mesh created by yams.
41 pos1=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
42 tip1=md.segments(pos1,1);
43 pos2=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));
44 tip2=md.segments(pos2,1);
45 if length(tip1)==2,
46 %swap.
47 temp=tip1; tip1=tip2; tip2=temp;
48 temp=pos1; pos1=pos2; pos2=temp;
49 pos=pos1;
50 else
51 pos=pos1;
52 end
53
54 pos_record=[pos_record; pos];
55 riftsegs=md.segments(pos,:);
56 while 1,
57 A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
58 %find other segment that holds B.
59 pos=find(md.segments(:,1)==B);
60 pos_record=[pos_record; pos];
61 riftsegs=[riftsegs; md.segments(pos,:)];
62 if ((riftsegs(end,2)==tip2(1)) | (riftsegs(end,2)==tip2(2))),
63 %figure out which tip we reached
64 if riftsegs(end,2)==tip2(1), index=2; else index=1; end
65 break;
66 end
67 end
68
69 %ok, now, we start from the other tip2, towards tip1
70 pos=pos2(index);
71 pos_record=[pos_record; pos];
72 riftsegs=[riftsegs; md.segments(pos,:)];
73 while 1,
74 A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
75 %find other segment that holds B.
76 pos=find(md.segments(:,1)==B);
77 pos_record=[pos_record; pos];
78 riftsegs=[riftsegs; md.segments(pos,:)];
79 if riftsegs(end,2)==tip1,
80 break;
81 end
82 end
83 md.rifts.riftstruct(i).segments=riftsegs;
84 md.rifts.riftstruct(i).tips=[tip1 tip2(1) tip2(2)];
85
86 end
87 end
88 end
89 %take out rift segments from segments
90 md.segments(pos_record,:)=[];
Note: See TracBrowser for help on using the repository browser.