- Timestamp:
- Apr 20, 2005, 5:41:03 PM (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution-parallel/general_mesh.py
r1178 r1237 6 6 7 7 A triangular element is defined in terms of three vertex ids, 8 ordered counter clock-wise, 8 ordered counter clock-wise, 9 9 each corresponding to a given coordinate set. 10 10 Vertices from different elements can point to the same … … 12 12 13 13 Coordinate sets are implemented as an N x 2 Numeric array containing 14 x and y coordinates. 15 14 x and y coordinates. 15 16 16 17 17 To instantiate: … … 33 33 c = [2.0,0.0] 34 34 e = [2.0, 2.0] 35 35 36 36 points = [a, b, c, e] 37 triangles = [ [1,0,2], [1,2,3] ] #bac, bce 37 triangles = [ [1,0,2], [1,2,3] ] #bac, bce 38 38 mesh = Mesh(points, triangles) 39 39 40 40 #creates two triangles: bac and bce 41 41 42 42 Other: 43 43 44 44 In addition mesh computes an Nx6 array called vertex_coordinates. 45 This structure is derived from coordinates and contains for each 45 This structure is derived from coordinates and contains for each 46 46 triangle the three x,y coordinates at the vertices. 47 47 48 48 49 49 This is a cut down version of mesh from pyvolution\mesh.py … … 58 58 59 59 origin is a 3-tuple consisting of UTM zone, easting and northing. 60 If specified coordinates are assumed to be relative to this origin. 60 If specified coordinates are assumed to be relative to this origin. 61 61 """ 62 62 … … 67 67 68 68 if geo_reference is None: 69 self.geo_reference = Geo_reference(DEFAULT_ZONE,0.0,0.0) 69 self.geo_reference = Geo_reference(DEFAULT_ZONE,0.0,0.0) 70 70 else: 71 71 self.geo_reference = geo_reference … … 76 76 77 77 msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples' 78 assert len(self.coordinates.shape) == 2, msg 78 assert len(self.coordinates.shape) == 2, msg 79 79 80 80 msg = 'Vertex indices reference non-existing coordinate sets' … … 84 84 #Register number of elements (N) 85 85 self.number_of_elements = N = self.triangles.shape[0] 86 87 #Allocate space for geometric quantities 86 87 #Allocate space for geometric quantities 88 88 self.normals = zeros((N, 6), Float) 89 89 self.areas = zeros(N, Float) … … 92 92 #Get x,y coordinates for all triangles and store 93 93 self.vertex_coordinates = V = self.compute_vertex_coordinates() 94 94 95 95 #Initialise each triangle 96 96 for i in range(N): 97 97 #if i % (N/10) == 0: print '(%d/%d)' %(i, N) 98 98 99 99 x0 = V[i, 0]; y0 = V[i, 1] 100 100 x1 = V[i, 2]; y1 = V[i, 3] 101 x2 = V[i, 4]; y2 = V[i, 5] 101 x2 = V[i, 4]; y2 = V[i, 5] 102 102 103 103 #Area … … 107 107 msg += ' is degenerate: area == %f' %self.areas[i] 108 108 assert self.areas[i] > 0.0, msg 109 109 110 110 111 111 #Normals … … 120 120 121 121 n0 = array([x2 - x1, y2 - y1]) 122 122 l0 = sqrt(sum(n0**2)) 123 123 124 124 n1 = array([x0 - x2, y0 - y2]) … … 137 137 n1[1], -n1[0], 138 138 n2[1], -n2[0]] 139 139 140 140 #Edgelengths 141 141 self.edgelengths[i, :] = [l0, l1, l2] 142 142 143 143 #Build vertex list 144 self.build_vertexlist() 145 144 self.build_vertexlist() 145 146 146 def __len__(self): 147 147 return self.number_of_elements … … 154 154 """Return all normal vectors. 155 155 Return normal vectors for all triangles as an Nx6 array 156 (ordered as x0, y0, x1, y1, x2, y2 for each triangle) 156 (ordered as x0, y0, x1, y1, x2, y2 for each triangle) 157 157 """ 158 158 return self.normals … … 163 163 Return value is the numeric array slice [x, y] 164 164 """ 165 return self.normals[i, 2*j:2*j+2] 166 167 168 165 return self.normals[i, 2*j:2*j+2] 166 167 168 169 169 def get_vertex_coordinates(self): 170 170 """Return all vertex coordinates. 171 171 Return all vertex coordinates for all triangles as an Nx6 array 172 (ordered as x0, y0, x1, y1, x2, y2 for each triangle) 172 (ordered as x0, y0, x1, y1, x2, y2 for each triangle) 173 173 """ 174 174 return self.vertex_coordinates … … 179 179 Return value is the numeric array slice [x, y] 180 180 """ 181 return self.vertex_coordinates[i, 2*j:2*j+2] 181 return self.vertex_coordinates[i, 2*j:2*j+2] 182 182 183 183 … … 189 189 #FIXME: Perhaps they should be ordered as in obj files?? 190 190 #See quantity.get_vertex_values 191 191 192 192 from Numeric import zeros, Float 193 193 194 194 N = self.number_of_elements 195 vertex_coordinates = zeros((N, 6), Float) 195 vertex_coordinates = zeros((N, 6), Float) 196 196 197 197 for i in range(N): … … 203 203 204 204 return vertex_coordinates 205 205 206 206 def get_vertices(self, indexes=None): 207 207 """Get connectivity 208 208 indexes is the set of element ids of interest 209 209 """ 210 210 211 211 from Numeric import take 212 212 213 213 if (indexes == None): 214 214 indexes = range(len(self)) #len(self)=number of elements 215 215 216 216 return take(self.triangles, indexes) 217 217 … … 220 220 unique_verts = {} 221 221 for triangle in triangles: 222 unique_verts[triangle[0]] = 0 223 unique_verts[triangle[1]] = 0 222 unique_verts[triangle[0]] = 0 223 unique_verts[triangle[1]] = 0 224 224 unique_verts[triangle[2]] = 0 225 return unique_verts.keys() 226 225 return unique_verts.keys() 226 227 227 def build_vertexlist(self): 228 228 """Build vertexlist index by vertex ids and for each entry (point id) … … 231 231 232 232 Preconditions: 233 self.coordinates and self.triangles are defined 234 233 self.coordinates and self.triangles are defined 234 235 235 Postcondition: 236 236 self.vertexlist is built … … 242 242 a = self.triangles[i, 0] 243 243 b = self.triangles[i, 1] 244 c = self.triangles[i, 2] 244 c = self.triangles[i, 2] 245 245 246 246 #Register the vertices v as lists of … … 251 251 vertexlist[v] = [] 252 252 253 vertexlist[v].append( (i, vertex_id) ) 253 vertexlist[v].append( (i, vertex_id) ) 254 254 255 255 self.vertexlist = vertexlist 256 256 257 257 258 258 def get_extent(self): … … 264 264 X = C[:,0:6:2].copy() 265 265 Y = C[:,1:6:2].copy() 266 266 267 267 xmin = min(X.flat) 268 268 xmax = max(X.flat) 269 269 ymin = min(Y.flat) 270 ymax = max(Y.flat) 271 270 ymax = max(Y.flat) 271 272 272 return xmin, xmax, ymin, ymax 273 273 … … 275 275 def get_area(self): 276 276 """Return total area of mesh 277 278 279 return sum(self.areas) 280 277 """ 278 279 return sum(self.areas) 280
Note: See TracChangeset
for help on using the changeset viewer.
