source: issm/trunk/src/m/exp/contourlevelzero.py@ 26744

Last change on this file since 26744 was 26744, checked in by Mathieu Morlighem, 3 years ago

merged trunk-jpl and trunk for revision 26742
File size: 8.9 KB
RevLine 
[26374]1import os.path
2import numpy as np
3from collections import OrderedDict
4
5def contourlevelzero(md,mask,level):
6 """CONTOURLEVELZERO - figure out the zero level (or offset thereof, specified by the level value)
7 of a vectorial mask, and vectorialize it into an exp or shp compatible structure.
8
9 Usage:
10 contours=contourlevelzero(md,mask,level)
11
12 See also: PLOT_CONTOUR
13 """
14
15 #process data
16 if md.mesh.dimension()==3:
17 x = md.mesh.x2d
18 y = md.mesh.y2d
19 z=md.mesh.z
20 index=md.mesh.elements2d-1
21 else:
22 x=md.mesh.x
23 y=md.mesh.y
24 index=md.mesh.elements-1
25 z=np.zeros((md.mesh.numberofvertices,1))
26
27 if len(mask)==0:
28 raise OSError("mask provided is empty")
29
30 if md.mesh.dimension()==3:
31 if len(mask)!=md.mesh.numberofvertices2d:
32 raise OSError("mask provided should be specified at the vertices of the mesh")
33 else:
34 if len(mask)!=md.mesh.numberofvertices:
35 raise OSError("mask provided should be specified at the vertices of the mesh")
36
37 #initialization of some variables
38 numberofelements=np.size(index,0)
39 elementslist=np.c_[0:numberofelements]
40 c=[]
41 h=[]
42
43 #get unique edges in mesh
44 #1: list of edges
45 edges=np.vstack((np.vstack((index[:,(0,1)],index[:,(1,2)])),index[:,(2,0)]))
46
47 #2: find unique edges
48 [edges,J]=np.unique(np.sort(edges,1),axis=0,return_inverse=True)
49 #3: unique edge numbers
50 vec=J
51 #4: unique edges numbers in each triangle (2 triangles sharing the same edge will have
52 # the same edge number)
53 edges_tria=np.hstack((np.hstack((vec[elementslist],vec[elementslist+numberofelements])),vec[elementslist+2*numberofelements]))
54
55 #segments [nodes1 nodes2]
56 Seg1=index[:,(0,1)]
57 Seg2=index[:,(1,2)]
58 Seg3=index[:,(2,0)]
59
60 #segment numbers [1;4;6;...]
61 Seg1_num=edges_tria[:,0]
62 Seg2_num=edges_tria[:,1]
63 Seg3_num=edges_tria[:,2]
64
65 #value of data on each tips of the segments
66 Data1=mask[Seg1]
67 Data2=mask[Seg2]
68 Data3=mask[Seg3]
69
70 #get the ranges for each segment
71 Range1=np.sort(Data1,1)
72 Range2=np.sort(Data2,1)
73 Range3=np.sort(Data3,1)
74
75 #find the segments that contain this value
76 pos1=(Range1[:,0]<level) & (Range1[:,1]>=level)
77 pos2=(Range2[:,0]<level) & (Range2[:,1]>=level)
78 pos3=(Range3[:,0]<level) & (Range3[:,1]>=level)
79
80 #get elements
81 poselem12=(pos1) & (pos2)
82 poselem13=(pos1) & (pos3)
83 poselem23=(pos2) & (pos3)
84 poselem=np.where((poselem12) | (poselem13) | (poselem23))
85 poselem=poselem[0]
86 numelems=len(poselem)
87
88 #if no element has been flagged, skip to the next level
89 if numelems==0:
90 raise Exception('contourlevelzero warning message: no elements found with corresponding level value in mask')
91 contours=[]
92 return contours
93
94 #go through the elements and build the coordinates for each segment (1 by element)
95 x1=np.zeros((numelems,1))
96 x2=np.zeros((numelems,1))
97 y1=np.zeros((numelems,1))
98 y2=np.zeros((numelems,1))
99 z1=np.zeros((numelems,1))
100 z2=np.zeros((numelems,1))
101
102 edge_l=np.zeros((numelems,2))
103
104 for j in range(0,numelems):
105
106 with np.errstate(divide='ignore', invalid='ignore'):
107 weight1=np.divide(level-Data1[poselem[j],0],Data1[poselem[j],1]-Data1[poselem[j],0])
108 weight2=np.divide(level-Data2[poselem[j],0],Data2[poselem[j],1]-Data2[poselem[j],0])
109 weight3=np.divide(level-Data3[poselem[j],0],Data3[poselem[j],1]-Data3[poselem[j],0])
110
111 if poselem12[poselem[j]]==True:
112
113 x1[j]=x[Seg1[poselem[j],0]]+weight1*[x[Seg1[poselem[j],1]]-x[Seg1[poselem[j],0]]]
114 x2[j]=x[Seg2[poselem[j],0]]+weight2*[x[Seg2[poselem[j],1]]-x[Seg2[poselem[j],0]]]
115 y1[j]=y[Seg1[poselem[j],0]]+weight1*[y[Seg1[poselem[j],1]]-y[Seg1[poselem[j],0]]]
116 y2[j]=y[Seg2[poselem[j],0]]+weight2*[y[Seg2[poselem[j],1]]-y[Seg2[poselem[j],0]]]
117 z1[j]=z[Seg1[poselem[j],0]]+weight1*[z[Seg1[poselem[j],1]]-z[Seg1[poselem[j],0]]]
118 z2[j]=z[Seg2[poselem[j],0]]+weight2*[z[Seg2[poselem[j],1]]-z[Seg2[poselem[j],0]]]
119
120 edge_l[j,0]=Seg1_num[poselem[j]]
121 edge_l[j,1]=Seg2_num[poselem[j]]
122 elif poselem13[poselem[j]]==True:
123
124 x1[j]=x[Seg1[poselem[j],0]]+weight1*[x[Seg1[poselem[j],1]]-x[Seg1[poselem[j],0]]]
125 x2[j]=x[Seg3[poselem[j],0]]+weight3*[x[Seg3[poselem[j],1]]-x[Seg3[poselem[j],0]]]
126 y1[j]=y[Seg1[poselem[j],0]]+weight1*[y[Seg1[poselem[j],1]]-y[Seg1[poselem[j],0]]]
127 y2[j]=y[Seg3[poselem[j],0]]+weight3*[y[Seg3[poselem[j],1]]-y[Seg3[poselem[j],0]]]
128 z1[j]=z[Seg1[poselem[j],0]]+weight1*[z[Seg1[poselem[j],1]]-z[Seg1[poselem[j],0]]]
129 z2[j]=z[Seg3[poselem[j],0]]+weight3*[z[Seg3[poselem[j],1]]-z[Seg3[poselem[j],0]]]
130
131 edge_l[j,0]=Seg1_num[poselem[j]]
132 edge_l[j,1]=Seg3_num[poselem[j]]
133 elif poselem23[poselem[j]]==True:
134
135 x1[j]=x[Seg2[poselem[j],0]]+weight2*[x[Seg2[poselem[j],1]]-x[Seg2[poselem[j],0]]]
136 x2[j]=x[Seg3[poselem[j],0]]+weight3*[x[Seg3[poselem[j],1]]-x[Seg3[poselem[j],0]]]
137 y1[j]=y[Seg2[poselem[j],0]]+weight2*[y[Seg2[poselem[j],1]]-y[Seg2[poselem[j],0]]]
138 y2[j]=y[Seg3[poselem[j],0]]+weight3*[y[Seg3[poselem[j],1]]-y[Seg3[poselem[j],0]]]
139 z1[j]=z[Seg2[poselem[j],0]]+weight2*[z[Seg2[poselem[j],1]]-z[Seg2[poselem[j],0]]]
140 z2[j]=z[Seg3[poselem[j],0]]+weight3*[z[Seg3[poselem[j],1]]-z[Seg3[poselem[j],0]]]
141
142 edge_l[j,0]=Seg2_num[poselem[j]]
143 edge_l[j,1]=Seg3_num[poselem[j]]
144
145 #else:
146 #it shoud not go here
147
148 #now that we have the segments, we must try to connect them...
149
150 #loop over the subcontours
151 contours=[]
152
153 while len(edge_l)>0:
154
155 #take the right edge of the second segment and connect it to the next segments if any
156 e1=edge_l[0,0]
157 e2=edge_l[0,1]
158 xc=np.vstack((x1[0],x2[0]))
159 yc=np.vstack((y1[0],y2[0]))
160 zc=np.vstack((z1[0],z2[0]))
161 #erase the lines corresponding to this edge
162 edge_l=np.delete(edge_l,0,axis=0)
163 x1=np.delete(x1,0,axis=0)
164 x2=np.delete(x2,0,axis=0)
165 y1=np.delete(y1,0,axis=0)
166 y2=np.delete(y2,0,axis=0)
167 z1=np.delete(z1,0,axis=0)
168 z2=np.delete(z2,0,axis=0)
169 pos1=np.where(edge_l==e1)
170
171 while len(pos1[0])>0:
172
173 if np.all(pos1[1]==0):
174 xc=np.vstack((x2[pos1[0]],xc))
175 yc=np.vstack((y2[pos1[0]],yc))
176 zc=np.vstack((z2[pos1[0]],zc))
177 #next edge:
178 e1=edge_l[pos1[0],1]
179 else:
180 xc=np.vstack((x1[pos1[0]],xc))
181 yc=np.vstack((y1[pos1[0]],yc))
182 zc=np.vstack((z1[pos1[0]],zc))
183 #next edge:
184 e1=edge_l[pos1[0],0]
185
186 #erase the lines of this
187 edge_l=np.delete(edge_l,pos1[0],axis=0)
188 x1=np.delete(x1,pos1[0],axis=0)
189 x2=np.delete(x2,pos1[0],axis=0)
190 y1=np.delete(y1,pos1[0],axis=0)
191 y2=np.delete(y2,pos1[0],axis=0)
192 z1=np.delete(z1,pos1[0],axis=0)
193 z2=np.delete(z2,pos1[0],axis=0)
194 #next connection
195 pos1=np.where(edge_l==e1)
196
197 #same thing the other way (to the right)
198 pos2=np.where(edge_l==e2)
199
200 while len(pos2[0])>0:
201
202 if np.all(pos2[1]==0):
203 xc=np.vstack((xc,x2[pos2[0]]))
204 yc=np.vstack((yc,y2[pos2[0]]))
205 zc=np.vstack((zc,z2[pos2[0]]))
206 #next edge:
207 e2=edge_l[pos2[0],1]
208 else:
209 xc=np.vstack((xc,x1[pos2[0]]))
210 yc=np.vstack((yc,y1[pos2[0]]))
211 zc=np.vstack((zc,z1[pos2[0]]))
212 #next edge:
213 e2=edge_l[pos2[0],0]
214
215 #erase the lines of this
216 edge_l=np.delete(edge_l,pos2[0],axis=0)
217 x1=np.delete(x1,pos2[0],axis=0)
218 x2=np.delete(x2,pos2[0],axis=0)
219 y1=np.delete(y1,pos2[0],axis=0)
220 y2=np.delete(y2,pos2[0],axis=0)
221 z1=np.delete(z1,pos2[0],axis=0)
222 z2=np.delete(z2,pos2[0],axis=0)
223 #next connection
224 pos2=np.where(edge_l==e2)
225
226 #save xc,yc contour:
227 newcontour = OrderedDict()
228 newcontour['nods'] = np.size(xc)
229 newcontour['density'] = 1
230 newcontour['closed'] = 0
231 newcontour['x'] = np.ma.filled(xc.astype(float), np.nan)
232 newcontour['y'] = np.ma.filled(yc.astype(float), np.nan)
233 newcontour['z'] = np.ma.filled(zc.astype(float), np.nan)
234 newcontour['name'] = ''
235 contours.append(newcontour)
236
237 return contours
Note: See TracBrowser for help on using the repository browser.