1 | """ |
---|
2 | This will produce a mesh which has vertex 10 that has no triangle. |
---|
3 | |
---|
4 | Why? |
---|
5 | |
---|
6 | Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou |
---|
7 | Geoscience Australia, 2004. |
---|
8 | """ |
---|
9 | |
---|
10 | |
---|
11 | import unittest |
---|
12 | from math import sqrt |
---|
13 | import time |
---|
14 | |
---|
15 | |
---|
16 | from pyvolution.least_squares import Interpolation |
---|
17 | from Numeric import allclose, array, transpose |
---|
18 | |
---|
19 | from coordinate_transforms.geo_reference import Geo_reference |
---|
20 | |
---|
21 | def distance(x, y): |
---|
22 | return sqrt( sum( (array(x)-array(y))**2 )) |
---|
23 | |
---|
24 | def linear_function(point): |
---|
25 | point = array(point) |
---|
26 | return point[:,0]+point[:,1] |
---|
27 | |
---|
28 | |
---|
29 | class Test_Least_Squares(unittest.TestCase): |
---|
30 | |
---|
31 | def setUp(self): |
---|
32 | pass |
---|
33 | |
---|
34 | def tearDown(self): |
---|
35 | pass |
---|
36 | |
---|
37 | def test_expand_search3(self): |
---|
38 | from random import seed, random |
---|
39 | from pmesh.mesh import Mesh |
---|
40 | |
---|
41 | num_of_points = 100 |
---|
42 | |
---|
43 | seed(2) |
---|
44 | m = Mesh() |
---|
45 | m.addUserVertex(0,0) |
---|
46 | m.addUserVertex(100,0) |
---|
47 | m.addUserVertex(0,100) |
---|
48 | m.addUserVertex(100,100) |
---|
49 | |
---|
50 | m.autoSegment(alpha = 100 ) |
---|
51 | |
---|
52 | dict = {} |
---|
53 | dict['points'] = [[10,10],[90,20]] |
---|
54 | dict['segments'] = [[0,1]] |
---|
55 | dict['segment_tags'] = ['wall1'] |
---|
56 | m.addVertsSegs(dict) |
---|
57 | |
---|
58 | dict = {} |
---|
59 | dict['points'] = [[10,90],[40,20]] |
---|
60 | dict['segments'] = [[0,1]] |
---|
61 | dict['segment_tags'] = ['wall2'] |
---|
62 | m.addVertsSegs(dict) |
---|
63 | |
---|
64 | dict = {} |
---|
65 | dict['points'] = [[20,90],[60,60]] |
---|
66 | dict['segments'] = [[0,1]] |
---|
67 | dict['segment_tags'] = ['wall3'] |
---|
68 | m.addVertsSegs(dict) |
---|
69 | |
---|
70 | dict = {} |
---|
71 | dict['points'] = [[60,20],[90,90]] |
---|
72 | dict['segments'] = [[0,1]] |
---|
73 | dict['segment_tags'] = ['wall4'] |
---|
74 | m.addVertsSegs(dict) |
---|
75 | |
---|
76 | m.addVertsSegs(dict) |
---|
77 | m.generateMesh(mode = "Q") |
---|
78 | m.export_mesh_file("aaaa.tsh") |
---|
79 | mesh_dict = m.Mesh2IOTriangulationDict() |
---|
80 | |
---|
81 | points = [] |
---|
82 | point_atts = [] |
---|
83 | for point in range(num_of_points): |
---|
84 | points.append([random()*100, random()*100]) |
---|
85 | point_atts.append(10.0) |
---|
86 | #m.autoSegment(alpha = 100 ) |
---|
87 | |
---|
88 | #print 'points', points |
---|
89 | t0 = time.time() |
---|
90 | interp = Interpolation(mesh_dict['vertices'], |
---|
91 | mesh_dict['triangles'], points, |
---|
92 | alpha=0.2, expand_search=True, |
---|
93 | verbose = True, |
---|
94 | max_points_per_cell = 4) |
---|
95 | calc = interp.fit_points(point_atts ) |
---|
96 | print "interp.expanded_quad_searches", interp.expanded_quad_searches |
---|
97 | #print "calc", calc |
---|
98 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
99 | |
---|
100 | #------------------------------------------------------------- |
---|
101 | if __name__ == "__main__": |
---|
102 | suite = unittest.makeSuite(Test_Least_Squares,'test') |
---|
103 | runner = unittest.TextTestRunner(verbosity=1) |
---|
104 | runner.run(suite) |
---|