Changeset 28028
- Timestamp:
- 12/19/23 13:19:39 (15 months ago)
- Location:
- issm/trunk-jpl/src/m/exp
- Files:
-
- 1 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/exp/contourlevelzero.py
r27950 r28028 1 import os.path2 import numpy as np3 from collections import OrderedDict4 5 6 1 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 2 raise Exception('this function has been renamed: A = isoline(md, mask)') -
issm/trunk-jpl/src/m/exp/isoline.m
r27961 r28028 166 166 edge_l(j,2)=Seg3_num(poselem(j)); 167 167 else 168 % it shoud not gohere168 %should never get here 169 169 end 170 170 end
Note:
See TracChangeset
for help on using the changeset viewer.