source: issm/oecreview/Archive/23390-24306/ISSM-23559-23560.diff@ 24307

Last change on this file since 24307 was 24307, checked in by Mathieu Morlighem, 5 years ago

NEW: adding Archive/23390-24306
File size: 6.1 KB
RevLine 
[24307]1Index: ../trunk-jpl/src/m/contrib/tsantos/mismip/if_position.m
2===================================================================
3--- ../trunk-jpl/src/m/contrib/tsantos/mismip/if_position.m (nonexistent)
4+++ ../trunk-jpl/src/m/contrib/tsantos/mismip/if_position.m (revision 23560)
5@@ -0,0 +1,199 @@
6+function [ifx ify] = if_position(md,step,level),
7+
8+ if(~isfield(md.results.TransientSolution,'MaskIceLevelset'))
9+ ifx=[];
10+ ify=[];
11+ return
12+ end
13+
14+ %initialization of some variables
15+ data = md.results.TransientSolution(step).MaskIceLevelset;
16+ if(isfield(md.results.TransientSolution,'MeshElements'))
17+ index = md.results.TransientSolution(step).MeshElements;
18+ x = md.results.TransientSolution(step).MeshX;
19+ y = md.results.TransientSolution(step).MeshY;
20+ else
21+ index = md.mesh.elements;
22+ x = md.mesh.x;
23+ y = md.mesh.y;
24+ end
25+ numberofelements = size(index,1);
26+ elementslist = 1:numberofelements;
27+ c = [];
28+ h = [];
29+
30+ %get unique edges in mesh
31+ %1: list of edges
32+ edges=[index(:,[1,2]); index(:,[2,3]); index(:,[3,1])];
33+ %2: find unique edges
34+ [edges,I,J]=unique(sort(edges,2),'rows');
35+ %3: unique edge numbers
36+ vec=J;
37+ %4: unique edges numbers in each triangle (2 triangles sharing the same edge will have
38+ % the same edge number)
39+ edges_tria=[vec(elementslist), vec(elementslist+numberofelements), vec(elementslist+2*numberofelements)];
40+
41+ %segments [nodes1 nodes2]
42+ Seg1=index(:,[1 2]);
43+ Seg2=index(:,[2 3]);
44+ Seg3=index(:,[3 1]);
45+
46+ %segment numbers [1;4;6;...]
47+ Seg1_num=edges_tria(:,1);
48+ Seg2_num=edges_tria(:,2);
49+ Seg3_num=edges_tria(:,3);
50+
51+ %value of data on each tips of the segments
52+ Data1=data(Seg1);
53+ Data2=data(Seg2);
54+ Data3=data(Seg3);
55+
56+ %get the ranges for each segment
57+ Range1=sort(Data1,2);
58+ Range2=sort(Data2,2);
59+ Range3=sort(Data3,2);
60+
61+ %find the segments that contain this value
62+ pos1=(Range1(:,1)<level & Range1(:,2)>level);
63+ pos2=(Range2(:,1)<level & Range2(:,2)>level);
64+ pos3=(Range3(:,1)<level & Range3(:,2)>level);
65+
66+ %get elements
67+ poselem12=(pos1 & pos2);
68+ poselem13=(pos1 & pos3);
69+ poselem23=(pos2 & pos3);
70+ poselem=find(poselem12 | poselem13 | poselem23);
71+ numelems=length(poselem);
72+
73+ %if no element has been flagged, skip to the next level
74+ if numelems==0,
75+ return,
76+ end
77+
78+ %go through the elements and build the coordinates for each segment (1 by element)
79+ x1=zeros(numelems,1);
80+ x2=zeros(numelems,1);
81+ y1=zeros(numelems,1);
82+ y2=zeros(numelems,1);
83+ edge_l=zeros(numelems,2);
84+
85+ for j=1:numelems,
86+
87+ weight1=(level-Data1(poselem(j),1))/(Data1(poselem(j),2)-Data1(poselem(j),1));
88+ weight2=(level-Data2(poselem(j),1))/(Data2(poselem(j),2)-Data2(poselem(j),1));
89+ weight3=(level-Data3(poselem(j),1))/(Data3(poselem(j),2)-Data3(poselem(j),1));
90+
91+ if poselem12(poselem(j));
92+
93+ x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1)));
94+ x2(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1)));
95+ y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
96+ y2(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
97+ edge_l(j,1)=Seg1_num(poselem(j));
98+ edge_l(j,2)=Seg2_num(poselem(j));
99+
100+ elseif poselem13(poselem(j)),
101+
102+ x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1)));
103+ x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1)));
104+ y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
105+ y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
106+ edge_l(j,1)=Seg1_num(poselem(j));
107+ edge_l(j,2)=Seg3_num(poselem(j));
108+
109+ elseif poselem23(poselem(j)),
110+
111+ x1(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1)));
112+ x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1)));
113+ y1(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
114+ y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
115+ edge_l(j,1)=Seg2_num(poselem(j));
116+ edge_l(j,2)=Seg3_num(poselem(j));
117+ else
118+ %it shoud not go here
119+ end
120+ end
121+
122+ %now that we have the segments, we must try to connect them...
123+
124+ %loop over the subcontours
125+ indice=0;
126+ while ~isempty(edge_l),
127+ indice=indice+1;
128+
129+ %take the right edge of the second segment and connect it to the next segments if any
130+ e1=edge_l(1,1); e2=edge_l(1,2);
131+ xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)];
132+
133+ %erase the lines corresponding to this edge
134+ edge_l(1,:)=[];
135+ x1(1)=[]; x2(1)=[];
136+ y1(1)=[]; y2(1)=[];
137+
138+ [ro1,co1]=find(edge_l==e1);
139+
140+ while ~isempty(ro1)
141+
142+ if co1==1,
143+ xc=[x2(ro1);xc]; yc=[y2(ro1);yc];
144+
145+ %next edge:
146+ e1=edge_l(ro1,2);
147+
148+ else
149+ xc=[x1(ro1);xc]; yc=[y1(ro1);yc];
150+
151+ %next edge:
152+ e1=edge_l(ro1,1);
153+ end
154+
155+ %erase the lines of this
156+ edge_l(ro1,:)=[];
157+ x1(ro1)=[]; x2(ro1)=[];
158+ y1(ro1)=[]; y2(ro1)=[];
159+
160+ %next connection
161+ [ro1,co1]=find(edge_l==e1);
162+ end
163+
164+ %same thing the other way (to the right)
165+ [ro2,co2]=find(edge_l==e2);
166+
167+ while ~isempty(ro2)
168+
169+ if co2==1,
170+ xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)];
171+
172+ %next edge:
173+ e2=edge_l(ro2,2);
174+ else
175+ xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)];
176+
177+ %next edge:
178+ e2=edge_l(ro2,1);
179+ end
180+
181+ %erase the lines of this
182+ edge_l(ro2,:)=[];
183+ x1(ro2)=[]; x2(ro2)=[];
184+ y1(ro2)=[]; y2(ro2)=[];
185+
186+ %next connection
187+ [ro2,co2]=find(edge_l==e2);
188+ end
189+
190+ % Update the CS data structure as per "contours.m"
191+ % so that clabel works
192+ c = horzcat(c,[level, xc'; length(xc), yc']);
193+ % y0=find(c(2,:)==0);
194+ % y50=find(c(2,:)==50000);
195+ % gl0(glstep)=c(1,y0);
196+ % gl50(glstep)=c(1,y50);
197+ ifx=c(1,1:end)';
198+ ify=c(2,1:end)';
199+ pos=find(ifx~=0);
200+ ifx=ifx(pos);
201+ ify=ify(pos);
202+
203+end
204+ %min(c(1,2:end))
205
206Property changes on: ../trunk-jpl/src/m/contrib/tsantos/mismip/if_position.m
207___________________________________________________________________
208Added: svn:executable
209## -0,0 +1 ##
210+*
211\ No newline at end of property
Note: See TracBrowser for help on using the repository browser.