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

Last change on this file was 28013, checked in by Mathieu Morlighem, 16 months ago

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