Changeset 3954
- Timestamp:
- Nov 9, 2006, 10:41:21 AM (18 years ago)
- Location:
- anuga_core/source
- Files:
-
- 6 edited
-
anuga/abstract_2d_finite_volumes/general_mesh.py (modified) (7 diffs)
-
anuga/abstract_2d_finite_volumes/quantity.py (modified) (5 diffs)
-
anuga/abstract_2d_finite_volumes/test_general_mesh.py (modified) (2 diffs)
-
anuga/shallow_water/data_manager.py (modified) (3 diffs)
-
anuga/shallow_water/shallow_water_domain.py (modified) (1 diff)
-
anuga_parallel/print_stats.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py
r3945 r3954 240 240 241 241 def get_nodes(self, absolute=False): 242 """Return all node coordinates ordered in an Nx2 array. 242 """Return all nodes in mesh. 243 244 The nodes are ordered in an Nx2 array where N is the number of nodes. 243 245 This is the same format they were provided in the constructor 244 246 i.e. without any duplication. … … 260 262 261 263 262 def get_vertex_coordinates(self, unique=False, absolute=False): 263 """Return all vertex coordinates. 264 Return all vertex coordinates for all triangles as a 3*N x 2 array 265 where the jth vertex of the ith triangle is located in row 3*i+j. 266 267 Boolean keyword unique will cause the points to be returned as 268 they were provided in the constructor i.e. without any duplication 269 in an N x 2 array. 270 271 """ 272 273 if unique is True: 274 return self.get_nodes(absolute) 275 276 277 V = self.vertex_coordinates 264 def get_vertex_coordinates(self, absolute=False): 265 """Return vertex coordinates for all triangles. 266 267 Return all vertex coordinates for all triangles as a 3*M x 2 array 268 where the jth vertex of the ith triangle is located in row 3*i+j and 269 M the number of triangles in the mesh. 270 271 Boolean keyword argument absolute determines whether coordinates 272 are to be made absolute by taking georeference into account 273 Default is False as many parts of ANUGA expects relative coordinates. 274 """ 275 276 277 278 V = self.vertex_coordinates #[:3*M,:] 278 279 if absolute is True: 279 280 if not self.geo_reference.is_absolute(): … … 294 295 295 296 def compute_vertex_coordinates(self): 296 """Return all vertex coordinates for all triangles as a 3* Nx 2 array297 """Return all vertex coordinates for all triangles as a 3*M x 2 array 297 298 where the jth vertex of the ith triangle is located in row 3*i+j. 298 299 … … 301 302 """ 302 303 303 N= self.number_of_triangles304 vertex_coordinates = zeros((3* N, 2), Float)305 306 for i in range( N):304 M = self.number_of_triangles 305 vertex_coordinates = zeros((3*M, 2), Float) 306 307 for i in range(M): 307 308 for j in range(3): 308 k = self.triangles[i,j] #Index of vertex 0 309 v_k = self.nodes[k] 310 311 vertex_coordinates[3*i+j,:] = v_k 312 309 k = self.triangles[i,j] # Index of vertex j in triangle i 310 vertex_coordinates[3*i+j,:] = self.nodes[k] 313 311 314 312 return vertex_coordinates 315 313 316 314 317 def get_vertices(self, indices=None): 318 """Get connectivity 319 indices is the set of element ids of interest 320 """ 321 322 N = self.number_of_full_triangles 315 316 def get_triangles(self, indices=None): 317 """Get mesh triangles. 318 319 Return Mx3 integer array where M is the number of triangles. 320 Each row corresponds to one triangle and the three entries are 321 indices into the mesh nodes which can be obtained using the method 322 get_nodes() 323 324 Optional argument, indices is the set of triangle ids of interest. 325 """ 326 327 M = self.number_of_full_triangles 323 328 324 329 if indices is None: 325 #indices = range(len(self)) #len(self)=number of elements 326 indices = range(N) 327 328 return take(self.triangles, indices) 330 indices = range(M) 331 332 return take(self.triangles, indices) 329 333 330 334 331 #FIXME - merge these two (get_vertices and get_triangles) 332 def get_triangles(self, obj=False): 333 """Get connetivity 334 Return triangles (triplets of indices into point coordinates) 335 336 If obj is True return structure commensurate with replicated 337 points, allowing for discontinuities 338 (FIXME: Need good name for this concept) 339 """ 340 341 if obj is True: 342 m = len(self) #Number of triangles 343 M = 3*m #Total number of unique vertices 344 T = reshape(array(range(M)).astype(Int), (m,3)) 345 else: 346 T = self.triangles 335 336 def get_disconnected_triangles(self): 337 """Get mesh based on nodes obtained from get_vertex_coordinates. 338 339 Return array Mx3 array of integers where each row corresponds to 340 a triangle. A triangle is a triplet of indices into 341 point coordinates obtained from get_vertex_coordinates and each 342 index appears only once 343 344 This provides a mesh where no triangles share nodes 345 (hence the name disconnected triangles) and different 346 nodes may have the same coordinates. 347 348 This version of the mesh is useful for storing meshes with 349 discontinuities at each node and is e.g. used for storing 350 data in sww files. 351 352 The triangles created will have the format 353 354 [[0,1,2], 355 [3,4,5], 356 [6,7,8], 357 ... 358 [3*M-3 3*M-2 3*M-1]] 359 """ 360 361 M = len(self) # Number of triangles 362 K = 3*M # Total number of unique vertices 363 T = reshape(array(range(K)).astype(Int), (M,3)) 347 364 348 365 return T … … 351 368 352 369 def get_unique_vertices(self, indices=None): 353 triangles = self.get_vertices(indices=indices) 370 """FIXME(Ole): This function needs a docstring 371 """ 372 triangles = self.get_triangles(indices=indices) 354 373 unique_verts = {} 355 374 for triangle in triangles: … … 361 380 362 381 def build_vertexlist(self): 363 """Build vertexlist index by vertex ids and for each entry (point id)382 """Build vertexlist indexed by vertex ids and for each entry (point id) 364 383 build a list of (triangles, vertex_id) pairs that use the point 365 384 as vertex. 366 385 386 The vertex list will have length N, where N is the number of nodes 387 in the mesh. 388 367 389 Preconditions: 368 390 self.nodes and self.triangles are defined … … 372 394 """ 373 395 374 vertexlist = [None]* len(self.nodes)396 vertexlist = [None]*self.number_of_nodes 375 397 for i in range(self.number_of_triangles): 376 398 -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
r3945 r3954 15 15 """ 16 16 17 from Numeric import array, zeros, Float, less, concatenate, NewAxis, argmax, allclose 17 from Numeric import array, zeros, Float, less, concatenate, NewAxis,\ 18 argmax, allclose 18 19 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar 19 20 … … 1026 1027 1027 1028 1028 #Method for outputting model results 1029 #FIXME: Split up into geometric and numeric stuff. 1030 #FIXME: Geometric (X,Y,V) should live in mesh.py 1031 #FIXME: STill remember to move XY to mesh 1029 # Methods for outputting model results 1032 1030 def get_vertex_values(self, 1033 1031 xy=True, 1034 smooth =None,1035 precision =None,1036 reduction =None):1037 """Return vertex values like an OBJ format 1032 smooth=None, 1033 precision=None, 1034 reduction=None): 1035 """Return vertex values like an OBJ format i.e. one value per node. 1038 1036 1039 1037 The vertex values are returned as one sequence in the 1D float array A. … … 1041 1039 1042 1040 The connectivity is represented as an integer array, V, of dimension 1043 M x 3, where M is the number of volumes. Each row has three indices 1044 into the X, Y, A arrays defining the triangle. 1041 Mx3, where M is the number of triangles. Each row has three indices 1042 defining the triangle and they correspond to elements in the arrays 1043 X, Y and A. 1045 1044 1046 1045 if smooth is True, vertex values corresponding to one common 1047 1046 coordinate set will be smoothed according to the given 1048 1047 reduction operator. In this case vertex coordinates will be 1049 de-duplicated. 1048 de-duplicated corresponding to the original nodes as obtained from 1049 the method general_mesh.get_nodes() 1050 1050 1051 1051 If no smoothings is required, vertex coordinates and values will 1052 1052 be aggregated as a concatenation of values at 1053 vertices 0, vertices 1 and vertices 2 1053 vertices 0, vertices 1 and vertices 2. This corresponds to 1054 the node coordinates obtained from the method 1055 general_mesh.get_vertex_coordinates() 1054 1056 1055 1057 … … 1066 1068 1067 1069 if smooth is None: 1070 # Take default from domain 1068 1071 smooth = self.domain.smooth 1069 1072 1070 1073 if precision is None: 1071 1074 precision = Float 1072 1073 #Create connectivity 1074 1075 if smooth == True: 1075 1076 1077 if smooth is True: 1078 # Ensure continuous vertex values by combining 1079 # values at each node using the reduction operator 1076 1080 1077 1081 if reduction is None: 1082 # Take default from domain 1078 1083 reduction = self.domain.reduction 1079 1084 1080 V = self.domain.get_ vertices()1081 N = len(self.domain.vertexlist)1085 V = self.domain.get_triangles() 1086 N = self.domain.number_of_full_nodes # Ignore ghost nodes if any 1082 1087 A = zeros(N, precision) 1083 1084 #Smoothing loop 1088 points = self.domain.get_nodes() 1089 1090 # Reduction loop 1085 1091 for k in range(N): 1086 1092 L = self.domain.vertexlist[k] 1087 1093 1088 #Go through all triangle, vertex pairs 1089 #contributing to vertex k and register vertex value 1090 1091 if L is None: continue #In case there are unused points 1092 1094 # Go through all triangle, vertex pairs 1095 # contributing to vertex k and register vertex value 1096 if L is None: continue # In case there are unused points 1093 1097 contributions = [] 1094 1098 for volume_id, vertex_id in L: … … 1098 1102 A[k] = reduction(contributions) 1099 1103 1100 1101 if xy is True: 1102 X = self.domain.get_nodes()[:,0].astype(precision) 1103 Y = self.domain.get_nodes()[:,1].astype(precision) 1104 1105 return X, Y, A, V 1106 else: 1107 return A, V 1108 else: 1109 #Don't smooth 1110 #obj machinery moved to general_mesh 1111 1112 # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]] 1113 # These vert_id's will relate to the verts created below 1114 #m = len(self.domain) #Number of volumes 1115 #M = 3*m #Total number of unique vertices 1116 #V = reshape(array(range(M)).astype(Int), (m,3)) 1117 1118 V = self.domain.get_triangles(obj=True) 1119 #FIXME use get_vertices, when ready 1120 1121 A = self.vertex_values.flat 1122 1123 #Do vertex coordinates 1124 if xy is True: 1125 C = self.domain.get_vertex_coordinates() 1126 1127 X = C[:,0:6:2].copy() 1128 Y = C[:,1:6:2].copy() 1129 1130 return X.flat, Y.flat, A, V 1131 else: 1132 return A, V 1104 else: 1105 # Allow discontinuous vertex values 1106 V = self.domain.get_disconnected_triangles() 1107 points = self.domain.get_vertex_coordinates() 1108 A = self.vertex_values.flat.astype(precision) 1109 1110 1111 # Return 1112 if xy is True: 1113 X = points[:,0].astype(precision) 1114 Y = points[:,1].astype(precision) 1115 1116 return X, Y, A, V 1117 else: 1118 return A, V 1119 1133 1120 1134 1121 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py
r3945 r3954 25 25 26 26 #Create basic mesh 27 points, vertices, boundary= rectangular(1, 3)28 domain = General_mesh( points, vertices)27 nodes, triangles, _ = rectangular(1, 3) 28 domain = General_mesh(nodes, triangles) 29 29 30 30 31 assert allclose(domain.get_vertex_coordinates(unique=True), 32 domain.nodes) 31 assert allclose(domain.get_nodes(), nodes) 33 32 34 #assert allclose(domain.get_vertex_coordinates(), ...TODO 35 #assert allclose(domain.get_vertex_coordinates(absolute=True), ...TODO 33 34 M = domain.number_of_triangles 35 36 V = domain.get_vertex_coordinates() 37 assert V.shape[0] == 3*M 38 39 for i in range(M): 40 for j in range(3): 41 k = triangles[i,j] #Index of vertex j in triangle i 42 assert allclose(V[3*i+j,:], nodes[k]) 43 44 36 45 37 46 … … 44 53 45 54 #Create basic mesh 46 points, vertices, boundary= rectangular(1, 3)47 domain = General_mesh( points, vertices)55 nodes, triangles, _ = rectangular(1, 3) 56 domain = General_mesh(nodes, triangles) 48 57 49 58 value = [7] 50 #indexes = [1] #FIXME (Ole): Should this be used 51 assert domain.get_vertices() == domain.triangles 52 assert domain.get_vertices([0,4]) == [domain.triangles[0], 53 domain.triangles[4]] 59 assert allclose(domain.get_triangles(), triangles) 60 assert allclose(domain.get_triangles([0,4]), 61 [triangles[0], triangles[4]]) 54 62 55 63 def test_areas(self): -
anuga_core/source/anuga/shallow_water/data_manager.py
r3953 r3954 372 372 precision=self.precision) 373 373 374 375 376 377 x[:] = X.astype(self.precision) 378 379 y[:] = Y.astype(self.precision) 380 381 # FIXME (HACK) 382 if len(z) <> len(Z): 383 truncation = self.domain.number_of_full_nodes 384 Z = Z[:truncation] 385 print len(z), len(Z), truncation, len(self.domain.get_nodes()) 386 387 z[:] = Z.astype(self.precision) 374 volumes[:] = V.astype(volumes.typecode()) 375 x[:] = X 376 y[:] = Y 377 z[:] = Z 388 378 389 379 #FIXME: Backwards compatibility 390 380 z = fid.variables['z'] 391 z[:] = Z .astype(self.precision)381 z[:] = Z 392 382 ################################ 393 383 394 volumes[:] = V.astype(volumes.typecode()) 384 395 385 396 386 #Close … … 501 491 precision = self.precision) 502 492 503 # HACK504 truncation = self.domain.number_of_full_nodes505 # HACK506 if stage.shape[1] <> len(A):507 A = A[:truncation]508 509 493 510 494 #FIXME: Make this general (see below) … … 2976 2960 outfile.variables['x'][:] = x - xllcorner 2977 2961 outfile.variables['y'][:] = y - yllcorner 2978 outfile.variables['z'][:] = z #FIXME HACK 2962 outfile.variables['z'][:] = z #FIXME HACK for bacwards compat. 2979 2963 outfile.variables['elevation'][:] = z 2980 2964 outfile.variables['time'][:] = times #Store time relative 2981 outfile.variables['volumes'][:] = volumes.astype(Int32) # OnOpteron 642965 outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64 2982 2966 2983 2967 -
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r3928 r3954 191 191 value per vertex using self.reduction. 192 192 """ 193 194 # FIXME (Ole): how about using the word continuous vertex values? 193 195 self.smooth = not flag 194 196 -
anuga_core/source/anuga_parallel/print_stats.py
r3514 r3954 46 46 def build_full_flag(domain, ghost_recv_dict): 47 47 48 tri_full_flag = ones(len(domain.get_ vertices()), Int8)48 tri_full_flag = ones(len(domain.get_triangles()), Int8) 49 49 for i in ghost_recv_dict.keys(): 50 50 for id in ghost_recv_dict[i][0]:
Note: See TracChangeset
for help on using the changeset viewer.
