source: trunk/anuga_work/development/analytical_solutions/fit_example.py @ 7924

"""Example where a surface is fitted to one set of points.
Then that surface is interpolated using another set of points.
"""

import sys
from os import sep
sys.path.append('..'+sep+'pyvolution')
from Numeric import array, allclose
from least_squares import Interpolation


#Simple function to provide example z values (z = x+y)
def linear_function(point):
    point = array(point)
    return point[:,0]+point[:,1]



#Setup mesh used to represent fitted function
a = [0.0, 0.0]
b = [0.0, 2.0]
c = [2.0, 0.0]
d = [0.0, 4.0]
e = [2.0, 2.0]
f = [4.0, 0.0]

points = [a, b, c, d, e, f]
#bac, bce, ecf, dbe, daf, dae 
triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]

#Datapoints to fit from
data_points1 = [[ 0.66666667, 0.66666667],
                [ 1.33333333, 1.33333333],
                [ 2.66666667, 0.66666667],
                [ 0.66666667, 2.66666667],
                [ 0.0, 1.0],
                [ 0.0, 3.0],
                [ 1.0, 0.0],
                [ 1.0, 1.0],
                [ 1.0, 2.0],
                [ 1.0, 3.0],                                          
                [ 2.0, 1.0],
                [ 3.0, 0.0],
                [ 3.0, 1.0]]

z = linear_function(data_points1) #Z-values


#Fit surface to mesh
interp = Interpolation(points, triangles, data_points1, alpha=0.0)
f = interp.fit(z)                 #Fitted values at vertices



#New datapoints where interpolated values are sought
data_points2 = [[ 0.0, 0.0],
                [ 0.5, 0.5],
                [ 1.0, 0.5],
                [ 2.8, 1.2]]                                          
        

#Build new A matrix based on new points
interp.build_interpolation_matrix_A(data_points2)

#Interpolate using fitted surface
z1 = interp.interpolate(f)

#Check result
answer = linear_function(data_points2)

print 'Result from fit:', z1
print 'Expected result:', answer

assert allclose(z1, answer)
print 'OK'