1 | import os.path
|
---|
2 | import numpy as np
|
---|
3 | from collections import OrderedDict
|
---|
4 |
|
---|
5 |
|
---|
6 | def 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
|
---|