- 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/test_neighbour_mesh.py
r6003 r6145 13 13 from mesh_factory import rectangular 14 14 from anuga.config import epsilon 15 from Numeric import allclose, array, Int16 15 17 16 from anuga.coordinate_transforms.geo_reference import Geo_reference … … 19 18 from anuga.utilities.numerical_tools import ensure_numeric 20 19 20 import Numeric as num 21 22 21 23 def distance(x, y): 22 return sqrt( sum( ( array(x)-array(y))**2 ))24 return sqrt( sum( (num.array(x)-num.array(y))**2 )) 23 25 24 26 class Test_Mesh(unittest.TestCase): … … 63 65 #Normals 64 66 normals = mesh.get_normals() 65 assert allclose(normals[0, 0:2], [3.0/5, 4.0/5])66 assert allclose(normals[0, 2:4], [-1.0, 0.0])67 assert allclose(normals[0, 4:6], [0.0, -1.0])68 69 assert allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])70 assert allclose(mesh.get_normal(0,1), [-1.0, 0.0])71 assert allclose(mesh.get_normal(0,2), [0.0, -1.0])67 assert num.allclose(normals[0, 0:2], [3.0/5, 4.0/5]) 68 assert num.allclose(normals[0, 2:4], [-1.0, 0.0]) 69 assert num.allclose(normals[0, 4:6], [0.0, -1.0]) 70 71 assert num.allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5]) 72 assert num.allclose(mesh.get_normal(0,1), [-1.0, 0.0]) 73 assert num.allclose(mesh.get_normal(0,2), [0.0, -1.0]) 72 74 73 75 #Edge lengths 74 assert allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])76 assert num.allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0]) 75 77 76 78 … … 81 83 82 84 V = mesh.get_vertex_coordinates() 83 assert allclose(V, [ [0.0, 0.0],85 assert num.allclose(V, [ [0.0, 0.0], 84 86 [4.0, 0.0], 85 87 [0.0, 3.0] ]) 86 88 87 89 V0 = mesh.get_vertex_coordinate(0, 0) 88 assert allclose(V0, [0.0, 0.0])90 assert num.allclose(V0, [0.0, 0.0]) 89 91 90 92 V1 = mesh.get_vertex_coordinate(0, 1) 91 assert allclose(V1, [4.0, 0.0])93 assert num.allclose(V1, [4.0, 0.0]) 92 94 93 95 V2 = mesh.get_vertex_coordinate(0, 2) 94 assert allclose(V2, [0.0, 3.0])96 assert num.allclose(V2, [0.0, 3.0]) 95 97 96 98 … … 237 239 238 240 mesh = Mesh(points, vertices,use_inscribed_circle=False) 239 assert allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work'241 assert num.allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work' 240 242 241 243 mesh = Mesh(points, vertices,use_inscribed_circle=True) 242 assert allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work'244 assert num.allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work' 243 245 244 246 def test_inscribed_circle_rightangle_triangle(self): … … 252 254 253 255 mesh = Mesh(points, vertices,use_inscribed_circle=False) 254 assert allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work'256 assert num.allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work' 255 257 256 258 mesh = Mesh(points, vertices,use_inscribed_circle=True) 257 assert allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work'259 assert num.allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work' 258 260 259 261 … … 269 271 assert mesh.areas[0] == 2.0 270 272 271 assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])273 assert num.allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3]) 272 274 273 275 … … 303 305 assert mesh.edgelengths[1,2] == sqrt(8.0) 304 306 305 assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])306 assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3])307 assert allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])308 assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])307 assert num.allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3]) 308 assert num.allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3]) 309 assert num.allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3]) 310 assert num.allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3]) 309 311 310 312 def test_mesh_and_neighbours(self): … … 698 700 """ 699 701 from mesh_factory import rectangular 700 from Numeric import zeros, Float701 702 702 703 #Create basic mesh … … 717 718 from mesh_factory import rectangular 718 719 #from mesh import Mesh 719 from Numeric import zeros, Float720 720 721 721 #Create basic mesh … … 727 727 728 728 assert len(P) == 8 729 assert allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],730 [1.0, 0.5], [1.0, 1.0], [0.5, 1.0],731 [0.0, 1.0], [0.0, 0.5]])729 assert num.allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0], 730 [1.0, 0.5], [1.0, 1.0], [0.5, 1.0], 731 [0.0, 1.0], [0.0, 0.5]]) 732 732 for p in points: 733 733 #print p, P … … 736 736 737 737 def test_boundary_polygon_II(self): 738 from Numeric import zeros, Float739 740 738 741 739 #Points … … 764 762 765 763 assert len(P) == 8 766 assert allclose(P, [a, d, f, i, h, e, c, b])764 assert num.allclose(P, [a, d, f, i, h, e, c, b]) 767 765 768 766 for p in points: … … 774 772 """Same as II but vertices ordered differently 775 773 """ 776 777 from Numeric import zeros, Float778 774 779 775 … … 804 800 805 801 assert len(P) == 8 806 assert allclose(P, [a, d, f, i, h, e, c, b])802 assert num.allclose(P, [a, d, f, i, h, e, c, b]) 807 803 808 804 for p in points: … … 815 811 is partitioned using pymetis. 816 812 """ 817 818 from Numeric import zeros, Float819 813 820 814 … … 847 841 848 842 # Note that point e appears twice! 849 assert allclose(P, [a, d, f, e, g, h, e, c, b])843 assert num.allclose(P, [a, d, f, e, g, h, e, c, b]) 850 844 851 845 for p in points: … … 864 858 """ 865 859 866 from Numeric import zeros, Float867 860 from mesh_factory import rectangular 868 861 … … 907 900 908 901 """ 909 from Numeric import zeros, Float910 911 902 912 903 #Points … … 955 946 956 947 assert len(P) == 8 957 assert allclose(P, [a, d, f, i, h, e, c, b])958 assert allclose(P, [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.5, 0.5), (1.0, 1.0), (0.5, 0.5), (0.0, 1.0), (0.0, 0.5)])948 assert num.allclose(P, [a, d, f, i, h, e, c, b]) 949 assert num.allclose(P, [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.5, 0.5), (1.0, 1.0), (0.5, 0.5), (0.0, 1.0), (0.0, 0.5)]) 959 950 960 951 … … 1101 1092 [ 35406.3359375 , 79332.9140625 ]] 1102 1093 1103 scaled_points = ensure_numeric(points, Int)/1000 # Simplify for ease of interpretation1094 scaled_points = ensure_numeric(points, num.Int)/1000 # Simplify for ease of interpretation 1104 1095 1105 1096 triangles = [[ 0, 1, 2], … … 1142 1133 assert is_inside_polygon(p, P) 1143 1134 1144 assert allclose(P, Pref)1135 assert num.allclose(P, Pref) 1145 1136 1146 1137 def test_lone_vertices(self): … … 1199 1190 boundary_polygon = mesh.get_boundary_polygon() 1200 1191 1201 assert allclose(absolute_points, boundary_polygon)1192 assert num.allclose(absolute_points, boundary_polygon) 1202 1193 1203 1194 def test_get_triangle_containing_point(self): … … 1246 1237 1247 1238 neighbours = mesh.get_triangle_neighbours(0) 1248 assert allclose(neighbours, [-1,1,-1])1239 assert num.allclose(neighbours, [-1,1,-1]) 1249 1240 neighbours = mesh.get_triangle_neighbours(-10) 1250 1241 assert neighbours == [] … … 1295 1286 for x in L: 1296 1287 if x.triangle_id % 2 == 0: 1297 assert allclose(x.length, ceiling-y_line)1288 assert num.allclose(x.length, ceiling-y_line) 1298 1289 else: 1299 assert allclose(x.length, y_line-floor)1290 assert num.allclose(x.length, y_line-floor) 1300 1291 1301 1292 1302 assert allclose(x.normal, [0,-1])1303 1304 assert allclose(x.segment[1][0], x.segment[0][0] + x.length)1305 assert allclose(x.segment[0][1], y_line)1306 assert allclose(x.segment[1][1], y_line)1293 assert num.allclose(x.normal, [0,-1]) 1294 1295 assert num.allclose(x.segment[1][0], x.segment[0][0] + x.length) 1296 assert num.allclose(x.segment[0][1], y_line) 1297 assert num.allclose(x.segment[1][1], y_line) 1307 1298 1308 1299 assert x.triangle_id in intersected_triangles … … 1311 1302 1312 1303 msg = 'Segments do not add up' 1313 assert allclose(total_length, 2), msg1304 assert num.allclose(total_length, 2), msg 1314 1305 1315 1306 … … 1346 1337 total_length = 0 1347 1338 for x in L: 1348 assert allclose(x.length, 1.0)1349 assert allclose(x.normal, [0,-1])1350 1351 assert allclose(x.segment[1][0], x.segment[0][0] + x.length)1352 assert allclose(x.segment[0][1], y_line)1353 assert allclose(x.segment[1][1], y_line)1339 assert num.allclose(x.length, 1.0) 1340 assert num.allclose(x.normal, [0,-1]) 1341 1342 assert num.allclose(x.segment[1][0], x.segment[0][0] + x.length) 1343 assert num.allclose(x.segment[0][1], y_line) 1344 assert num.allclose(x.segment[1][1], y_line) 1354 1345 1355 1346 … … 1360 1351 1361 1352 msg = 'Segments do not add up' 1362 assert allclose(total_length, 2), msg1353 assert num.allclose(total_length, 2), msg 1363 1354 1364 1355 … … 1400 1391 for x in L: 1401 1392 if x.triangle_id == 1: 1402 assert allclose(x.length, 1)1403 assert allclose(x.normal, [0, -1])1393 assert num.allclose(x.length, 1) 1394 assert num.allclose(x.normal, [0, -1]) 1404 1395 1405 1396 if x.triangle_id == 5: 1406 assert allclose(x.length, 0.5)1407 assert allclose(x.normal, [0, -1])1397 assert num.allclose(x.length, 0.5) 1398 assert num.allclose(x.normal, [0, -1]) 1408 1399 1409 1400 … … 1413 1404 1414 1405 msg = 'Segments do not add up' 1415 assert allclose(total_length, 1.5), msg1406 assert num.allclose(total_length, 1.5), msg 1416 1407 1417 1408 … … 1447 1438 total_length = 0 1448 1439 for i, x in enumerate(L): 1449 assert allclose(x.length, s2)1450 assert allclose(x.normal, [-s2, -s2])1451 assert allclose(sum(x.normal**2), 1)1440 assert num.allclose(x.length, s2) 1441 assert num.allclose(x.normal, [-s2, -s2]) 1442 assert num.allclose(sum(x.normal**2), 1) 1452 1443 1453 1444 assert x.triangle_id in intersected_triangles … … 1456 1447 1457 1448 msg = 'Segments do not add up' 1458 assert allclose(total_length, 4*s2), msg1449 assert num.allclose(total_length, 4*s2), msg 1459 1450 1460 1451 … … 1471 1462 total_length = 0 1472 1463 for i, x in enumerate(L): 1473 assert allclose(x.length, s2)1474 assert allclose(x.normal, [s2, s2])1475 assert allclose(sum(x.normal**2), 1)1464 assert num.allclose(x.length, s2) 1465 assert num.allclose(x.normal, [s2, s2]) 1466 assert num.allclose(sum(x.normal**2), 1) 1476 1467 1477 1468 assert x.triangle_id in intersected_triangles … … 1480 1471 1481 1472 msg = 'Segments do not add up' 1482 assert allclose(total_length, 4*s2), msg1473 assert num.allclose(total_length, 4*s2), msg 1483 1474 1484 1475 … … 1496 1487 total_length = 0 1497 1488 for i, x in enumerate(L): 1498 assert allclose(x.length, 2*s2)1499 assert allclose(x.normal, [-s2, s2])1500 assert allclose(sum(x.normal**2), 1)1489 assert num.allclose(x.length, 2*s2) 1490 assert num.allclose(x.normal, [-s2, s2]) 1491 assert num.allclose(sum(x.normal**2), 1) 1501 1492 1502 1493 assert x.triangle_id in intersected_triangles … … 1505 1496 1506 1497 msg = 'Segments do not add up' 1507 assert allclose(total_length, 4*s2), msg1498 assert num.allclose(total_length, 4*s2), msg 1508 1499 1509 1500 … … 1520 1511 total_length = 0 1521 1512 for i, x in enumerate(L): 1522 assert allclose(x.length, 2*s2)1523 assert allclose(x.normal, [s2, -s2])1524 assert allclose(sum(x.normal**2), 1)1513 assert num.allclose(x.length, 2*s2) 1514 assert num.allclose(x.normal, [s2, -s2]) 1515 assert num.allclose(sum(x.normal**2), 1) 1525 1516 1526 1517 assert x.triangle_id in intersected_triangles … … 1529 1520 1530 1521 msg = 'Segments do not add up' 1531 assert allclose(total_length, 4*s2), msg1522 assert num.allclose(total_length, 4*s2), msg 1532 1523 1533 1524 … … 1545 1536 total_length = 0 1546 1537 for i, x in enumerate(L): 1547 assert allclose(x.length, s2)1548 assert allclose(x.normal, [-s2, -s2])1549 assert allclose(sum(x.normal**2), 1)1538 assert num.allclose(x.length, s2) 1539 assert num.allclose(x.normal, [-s2, -s2]) 1540 assert num.allclose(sum(x.normal**2), 1) 1550 1541 1551 1542 assert x.triangle_id in intersected_triangles … … 1554 1545 1555 1546 msg = 'Segments do not add up' 1556 assert allclose(total_length, 2*s2), msg1547 assert num.allclose(total_length, 2*s2), msg 1557 1548 1558 1549 … … 1567 1558 total_length = 0 1568 1559 for i, x in enumerate(L): 1569 assert allclose(x.normal, [-s2, -s2])1570 assert allclose(sum(x.normal**2), 1)1560 assert num.allclose(x.normal, [-s2, -s2]) 1561 assert num.allclose(sum(x.normal**2), 1) 1571 1562 1572 1563 msg = 'Triangle %d' %x.triangle_id + ' is not in %s' %(intersected_triangles) … … 1603 1594 assert len(L) == 1 1604 1595 assert L[0].triangle_id == 3 1605 assert allclose(L[0].length, 0.5)1606 assert allclose(L[0].normal, [-1,0])1596 assert num.allclose(L[0].length, 0.5) 1597 assert num.allclose(L[0].normal, [-1,0]) 1607 1598 1608 1599 … … 1615 1606 assert len(L) == 1 1616 1607 assert L[0].triangle_id == 3 1617 assert allclose(L[0].length, 0.4)1618 assert allclose(L[0].normal, [-1,0])1608 assert num.allclose(L[0].length, 0.4) 1609 assert num.allclose(L[0].normal, [-1,0]) 1619 1610 1620 1611 intersected_triangles = [3] … … 1628 1619 assert len(L) == 1 1629 1620 assert L[0].triangle_id == 3 1630 assert allclose(L[0].length, 0.4)1631 assert allclose(L[0].normal, [1,0])1621 assert num.allclose(L[0].length, 0.4) 1622 assert num.allclose(L[0].normal, [1,0]) 1632 1623 1633 1624 … … 1663 1654 for x in L: 1664 1655 if x.triangle_id == 3: 1665 assert allclose(x.length, 0.5)1666 assert allclose(x.normal, [-1,0])1656 assert num.allclose(x.length, 0.5) 1657 assert num.allclose(x.normal, [-1,0]) 1667 1658 1668 1659 if x.triangle_id == 2: 1669 assert allclose(x.length, s2)1670 assert allclose(x.normal, [-s2,-s2])1660 assert num.allclose(x.length, s2) 1661 assert num.allclose(x.normal, [-s2,-s2]) 1671 1662 1672 1663 … … 1701 1692 for x in L: 1702 1693 if x.triangle_id == 3: 1703 assert allclose(x.length, 0.5)1704 assert allclose(x.normal, [-1,0])1694 assert num.allclose(x.length, 0.5) 1695 assert num.allclose(x.normal, [-1,0]) 1705 1696 1706 1697 if x.triangle_id == 2: 1707 1698 msg = str(x.length) 1708 assert allclose(x.length, s2), msg1709 assert allclose(x.normal, [-s2,-s2])1699 assert num.allclose(x.length, s2), msg 1700 assert num.allclose(x.normal, [-s2,-s2]) 1710 1701 1711 1702 if x.triangle_id == 5: 1712 segvec = array([line[2][0]-1,1713 line[2][1]-1])1703 segvec = num.array([line[2][0]-1, 1704 line[2][1]-1]) 1714 1705 msg = str(x.length) 1715 assert allclose(x.length, sqrt(sum(segvec**2))), msg1716 assert allclose(x.normal, [-s2,-s2])1706 assert num.allclose(x.length, sqrt(sum(segvec**2))), msg 1707 assert num.allclose(x.normal, [-s2,-s2]) 1717 1708 1718 1709 … … 1752 1743 for x in L: 1753 1744 if x.triangle_id == 3: 1754 assert allclose(x.length, 0.5)1755 assert allclose(x.normal, [-1,0])1745 assert num.allclose(x.length, 0.5) 1746 assert num.allclose(x.normal, [-1,0]) 1756 1747 1757 1748 if x.triangle_id == 2: 1758 1749 msg = str(x.length) 1759 assert allclose(x.length, s2), msg1760 assert allclose(x.normal, [-s2,-s2])1750 assert num.allclose(x.length, s2), msg 1751 assert num.allclose(x.normal, [-s2,-s2]) 1761 1752 1762 1753 if x.triangle_id == 5: 1763 1754 if x.segment == ((1.0, 1.0), (1.25, 0.75)): 1764 segvec = array([line[2][0]-1,1765 line[2][1]-1])1755 segvec = num.array([line[2][0]-1, 1756 line[2][1]-1]) 1766 1757 msg = str(x.length) 1767 assert allclose(x.length, sqrt(sum(segvec**2))), msg1768 assert allclose(x.normal, [-s2,-s2])1758 assert num.allclose(x.length, sqrt(sum(segvec**2))), msg 1759 assert num.allclose(x.normal, [-s2,-s2]) 1769 1760 elif x.segment == ((1.25, 0.75), (1.5, 1.0)): 1770 segvec = array([1.5-line[2][0],1771 1.0-line[2][1]])1761 segvec = num.array([1.5-line[2][0], 1762 1.0-line[2][1]]) 1772 1763 1773 assert allclose(x.length, sqrt(sum(segvec**2))), msg1774 assert allclose(x.normal, [s2,-s2])1764 assert num.allclose(x.length, sqrt(sum(segvec**2))), msg 1765 assert num.allclose(x.normal, [s2,-s2]) 1775 1766 else: 1776 1767 msg = 'Unknow segment: %s' %x.segment … … 1780 1771 1781 1772 if x.triangle_id == 6: 1782 assert allclose(x.normal, [s2,-s2])1783 assert allclose(x.segment, ((1.5, 1.0), (2, 1.5)))1773 assert num.allclose(x.normal, [s2,-s2]) 1774 assert num.allclose(x.segment, ((1.5, 1.0), (2, 1.5))) 1784 1775 1785 1776 … … 1791 1782 #xi1 = line[1][0] 1792 1783 #eta1 = line[1][1] 1793 #linevector = array([xi1-xi0, eta1-eta0])1784 #linevector = num.array([xi1-xi0, eta1-eta0]) 1794 1785 #linelength = sqrt(sum(linevector**2)) 1795 1786 # … … 1847 1838 ref_length = line[1][1] - line[0][1] 1848 1839 #print ref_length, total_length 1849 assert allclose(total_length, ref_length)1840 assert num.allclose(total_length, ref_length) 1850 1841 1851 1842
Note: See TracChangeset
for help on using the changeset viewer.