- Timestamp:
- Jan 13, 2009, 11:51:22 AM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
r5897 r6145 11 11 from anuga.caching import cache 12 12 from math import pi, sqrt 13 from Numeric import array, allclose 13 14 import Numeric as num 14 15 15 16 … … 81 82 82 83 83 from Numeric import array, zeros, Int, Float, maximum, sqrt, sum84 84 85 85 General_mesh.__init__(self, coordinates, triangles, … … 98 98 99 99 #Allocate space for geometric quantities 100 self.centroid_coordinates = zeros((N, 2),Float)101 102 self.radii = zeros(N,Float)103 104 self.neighbours = zeros((N, 3),Int)105 self.neighbour_edges = zeros((N, 3),Int)106 self.number_of_boundaries = zeros(N,Int)107 self.surrogate_neighbours = zeros((N, 3),Int)100 self.centroid_coordinates = num.zeros((N, 2), num.Float) 101 102 self.radii = num.zeros(N, num.Float) 103 104 self.neighbours = num.zeros((N, 3), num.Int) 105 self.neighbour_edges = num.zeros((N, 3), num.Int) 106 self.number_of_boundaries = num.zeros(N, num.Int) 107 self.surrogate_neighbours = num.zeros((N, 3), num.Int) 108 108 109 109 #Get x,y coordinates for all triangles and store … … 124 124 125 125 #Compute centroid 126 centroid = array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])126 centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3]) 127 127 self.centroid_coordinates[i] = centroid 128 128 … … 133 133 134 134 #Midpoints 135 m0 = array([(x1 + x2)/2, (y1 + y2)/2])136 m1 = array([(x0 + x2)/2, (y0 + y2)/2])137 m2 = array([(x1 + x0)/2, (y1 + y0)/2])135 m0 = num.array([(x1 + x2)/2, (y1 + y2)/2]) 136 m1 = num.array([(x0 + x2)/2, (y0 + y2)/2]) 137 m2 = num.array([(x1 + x0)/2, (y1 + y0)/2]) 138 138 139 139 #The radius is the distance from the centroid of 140 140 #a triangle to the midpoint of the side of the triangle 141 141 #closest to the centroid 142 d0 = sqrt(sum( (centroid-m0)**2 ))143 d1 = sqrt(sum( (centroid-m1)**2 ))144 d2 = sqrt(sum( (centroid-m2)**2 ))142 d0 = num.sqrt(num.sum( (centroid-m0)**2 )) 143 d1 = num.sqrt(num.sum( (centroid-m1)**2 )) 144 d2 = num.sqrt(num.sum( (centroid-m2)**2 )) 145 145 146 146 self.radii[i] = min(d0, d1, d2) … … 150 150 #of inscribed circle is computed 151 151 152 a = sqrt((x0-x1)**2+(y0-y1)**2)153 b = sqrt((x1-x2)**2+(y1-y2)**2)154 c = sqrt((x2-x0)**2+(y2-y0)**2)152 a = num.sqrt((x0-x1)**2+(y0-y1)**2) 153 b = num.sqrt((x1-x2)**2+(y1-y2)**2) 154 c = num.sqrt((x2-x0)**2+(y2-y0)**2) 155 155 156 156 self.radii[i]=2.0*self.areas[i]/(a+b+c) … … 212 212 x1 = V[i, 2]; y1 = V[i, 3] 213 213 x2 = V[i, 4]; y2 = V[i, 5] 214 a = sqrt((x0-x1)**2+(y0-y1)**2)215 b = sqrt((x1-x2)**2+(y1-y2)**2)216 c = sqrt((x2-x0)**2+(y2-y0)**2)214 a = num.sqrt((x0-x1)**2+(y0-y1)**2) 215 b = num.sqrt((x1-x2)**2+(y1-y2)**2) 216 c = num.sqrt((x2-x0)**2+(y2-y0)**2) 217 217 ratio = old_rad/self.radii[i] 218 218 max_ratio = ratio … … 224 224 x1 = V[i, 2]; y1 = V[i, 3] 225 225 x2 = V[i, 4]; y2 = V[i, 5] 226 a = sqrt((x0-x1)**2+(y0-y1)**2)227 b = sqrt((x1-x2)**2+(y1-y2)**2)228 c = sqrt((x2-x0)**2+(y2-y0)**2)226 a = num.sqrt((x0-x1)**2+(y0-y1)**2) 227 b = num.sqrt((x1-x2)**2+(y1-y2)**2) 228 c = num.sqrt((x2-x0)**2+(y2-y0)**2) 229 229 self.radii[i]=self.areas[i]/(2*(a+b+c))*safety_factor 230 230 ratio = old_rad/self.radii[i] … … 388 388 self.element_tag is defined 389 389 """ 390 from Numeric import array, Int391 390 392 391 if tagged_elements is None: … … 395 394 #Check that all keys in given boundary exist 396 395 for tag in tagged_elements.keys(): 397 tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)396 tagged_elements[tag] = num.array(tagged_elements[tag]).astype(num.Int) 398 397 399 398 msg = 'Not all elements exist. ' … … 463 462 """ 464 463 465 from Numeric import allclose, sqrt, array, minimum, maximum466 464 from anuga.utilities.numerical_tools import angle, ensure_numeric 467 465 … … 477 475 inverse_segments = {} 478 476 p0 = None 479 mindist = sqrt(sum((pmax-pmin)**2)) # Start value across entire mesh477 mindist = num.sqrt(num.sum((pmax-pmin)**2)) # Start value across entire mesh 480 478 for i, edge_id in self.boundary.keys(): 481 479 # Find vertex ids for boundary segment … … 491 489 # Note: Could be arbitrary, but nice to have 492 490 # a unique way of selecting 493 dist_A = sqrt(sum((A-pmin)**2))494 dist_B = sqrt(sum((B-pmin)**2))491 dist_A = num.sqrt(num.sum((A-pmin)**2)) 492 dist_B = num.sqrt(num.sum((B-pmin)**2)) 495 493 496 494 # Find lower leftmost point … … 593 591 # We have reached a point already visited. 594 592 595 if allclose(p1, polygon[0]):593 if num.allclose(p1, polygon[0]): 596 594 # If it is the initial point, the polygon is complete. 597 595 … … 629 627 from anuga.utilities.numerical_tools import anglediff 630 628 631 from Numeric import sort, allclose632 633 629 N = len(self) 634 630 … … 672 668 673 669 # Normalise 674 l_u = sqrt(u[0]*u[0] + u[1]*u[1])675 l_v = sqrt(v[0]*v[0] + v[1]*v[1])670 l_u = num.sqrt(u[0]*u[0] + u[1]*u[1]) 671 l_v = num.sqrt(v[0]*v[0] + v[1]*v[1]) 676 672 677 673 msg = 'Normal vector in triangle %d does not have unit length' %i 678 assert allclose(l_v, 1), msg674 assert num.allclose(l_v, 1), msg 679 675 680 676 x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product … … 730 726 731 727 V = self.vertex_value_indices[:] #Take a copy 732 V = sort(V)733 assert allclose(V, range(3*N))734 735 assert sum(self.number_of_triangles_per_node) ==\736 len(self.vertex_value_indices)728 V = num.sort(V) 729 assert num.allclose(V, range(3*N)) 730 731 assert num.sum(self.number_of_triangles_per_node) ==\ 732 len(self.vertex_value_indices) 737 733 738 734 # Check number of triangles per node … … 742 738 count[i] += 1 743 739 744 assert allclose(count, self.number_of_triangles_per_node)740 assert num.allclose(count, self.number_of_triangles_per_node) 745 741 746 742 … … 805 801 """ 806 802 807 from Numeric import arange808 803 from anuga.utilities.numerical_tools import histogram, create_bins 809 804 … … 1141 1136 1142 1137 # Distances from line origin to the two intersections 1143 z0 = array([x0 - xi0, y0 - eta0])1144 z1 = array([x1 - xi0, y1 - eta0])1145 d0 = sqrt(sum(z0**2))1146 d1 = sqrt(sum(z1**2))1138 z0 = num.array([x0 - xi0, y0 - eta0]) 1139 z1 = num.array([x1 - xi0, y1 - eta0]) 1140 d0 = num.sqrt(num.sum(z0**2)) 1141 d1 = num.sqrt(num.sum(z1**2)) 1147 1142 1148 1143 if d1 < d0: … … 1157 1152 # Normal direction: 1158 1153 # Right hand side relative to line direction 1159 vector = array([x1 - x0, y1 - y0]) # Segment vector1160 length = sqrt(sum(vector**2)) # Segment length1161 normal = array([vector[1], -vector[0]])/length1154 vector = num.array([x1 - x0, y1 - y0]) # Segment vector 1155 length = num.sqrt(num.sum(vector**2)) # Segment length 1156 normal = num.array([vector[1], -vector[0]])/length 1162 1157 1163 1158 … … 1239 1234 assert isinstance(segment, Triangle_intersection), msg 1240 1235 1241 midpoint = sum(array(segment.segment))/21236 midpoint = num.sum(num.array(segment.segment))/2 1242 1237 midpoints.append(midpoint) 1243 1238
Note: See TracChangeset
for help on using the changeset viewer.