source: inundation/ga/storm_surge/analytical solutions/fit_example.py @ 643

Last change on this file since 643 was 643, checked in by ole, 20 years ago
File size: 1.9 KB
Line 
1
2"""Example where a surface is fitted to one set of points.
3Then that surface is interpolated using another set of points.
4"""
5
6import sys
7from os import sep
8sys.path.append('..'+sep+'pyvolution')
9from Numeric import array, allclose
10from least_squares import Interpolation
11
12
13#Simple function to provide example z values (z = x+y)
14def linear_function(point):
15    point = array(point)
16    return point[:,0]+point[:,1]
17
18
19
20#Setup mesh used to represent fitted function
21a = [0.0, 0.0]
22b = [0.0, 2.0]
23c = [2.0, 0.0]
24d = [0.0, 4.0]
25e = [2.0, 2.0]
26f = [4.0, 0.0]
27
28points = [a, b, c, d, e, f]
29#bac, bce, ecf, dbe, daf, dae
30triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
31
32#Datapoints to fit from
33data_points1 = [[ 0.66666667, 0.66666667],
34                [ 1.33333333, 1.33333333],
35                [ 2.66666667, 0.66666667],
36                [ 0.66666667, 2.66666667],
37                [ 0.0, 1.0],
38                [ 0.0, 3.0],
39                [ 1.0, 0.0],
40                [ 1.0, 1.0],
41                [ 1.0, 2.0],
42                [ 1.0, 3.0],                                         
43                [ 2.0, 1.0],
44                [ 3.0, 0.0],
45                [ 3.0, 1.0]]
46
47z = linear_function(data_points1) #Z-values
48
49
50#Fit surface to mesh
51interp = Interpolation(points, triangles, data_points1, alpha=0.0)
52f = interp.fit(z)                 #Fitted values at vertices
53
54
55
56#New datapoints where interpolated values are sought
57data_points2 = [[ 0.0, 0.0],
58                [ 0.5, 0.5],
59                [ 1.0, 0.5],
60                [ 2.8, 1.2]]                                         
61       
62
63#Build new A matrix based on new points
64interp.build_interpolation_matrix_A(data_points2)
65
66#Interpolate using fitted surface
67z1 = interp.interpolate(f)
68
69#Check result
70answer = linear_function(data_points2)
71
72print 'Result from fit:', z1
73print 'Expected result:', answer
74
75assert allclose(z1, answer)
76print 'OK'
Note: See TracBrowser for help on using the repository browser.