Context Navigation

source: inundation/fit_interpolate/test_interpolate.py @ 3415

Last change on this file since 3415 was 3412, checked in by duncan, 19 years ago
typo
File size: 49.7 KB

Rev	Line
[2187]	1	#!/usr/bin/env python
	2
	3	#TEST
[2394]	4
	5	#import time, os
	6
	7
[2187]	8	import sys
[2394]	9	import os
[2187]	10	import unittest
	11	from math import sqrt
[2394]	12	import tempfile
[2187]	13
[2394]	14	from Scientific.IO.NetCDF import NetCDFFile
	15	from Numeric import allclose, array, transpose, zeros, Float
[2187]	16
[2394]	17
	18	# ANUGA code imports
[2187]	19	from interpolate import *
	20	from coordinate_transforms.geo_reference import Geo_reference
[2651]	21	from shallow_water import Domain, Transmissive_boundary
[2655]	22	from utilities.numerical_tools import mean, INF
[2394]	23	from data_manager import get_dataobject
[2577]	24	from geospatial_data.geospatial_data import Geospatial_data
[3411]	25	from pmesh.mesh import Mesh
[2187]	26
	27	def distance(x, y):
	28	return sqrt( sum( (array(x)-array(y))**2 ))
	29
	30	def linear_function(point):
	31	point = array(point)
	32	return point[:,0]+point[:,1]
	33
	34
	35	class Test_Interpolate(unittest.TestCase):
	36
	37	def setUp(self):
	38
[2394]	39	import time
	40	from mesh_factory import rectangular
	41
	42
	43	#Create basic mesh
	44	points, vertices, boundary = rectangular(2, 2)
	45
	46	#Create shallow water domain
	47	domain = Domain(points, vertices, boundary)
	48	domain.default_order=2
	49	domain.beta_h = 0
	50
	51
	52	#Set some field values
	53	domain.set_quantity('elevation', lambda x,y: -x)
	54	domain.set_quantity('friction', 0.03)
	55
	56
	57	######################
	58	# Boundary conditions
	59	B = Transmissive_boundary(domain)
	60	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
	61
	62
	63	######################
	64	#Initial condition - with jumps
	65
	66	bed = domain.quantities['elevation'].vertex_values
	67	stage = zeros(bed.shape, Float)
	68
	69	h = 0.3
	70	for i in range(stage.shape[0]):
	71	if i % 2 == 0:
	72	stage[i,:] = bed[i,:] + h
	73	else:
	74	stage[i,:] = bed[i,:]
	75
	76	domain.set_quantity('stage', stage)
	77
	78	domain.distribute_to_vertices_and_edges()
	79
	80
	81	self.domain = domain
	82
	83	C = domain.get_vertex_coordinates()
	84	self.X = C[:,0:6:2].copy()
	85	self.Y = C[:,1:6:2].copy()
	86
	87	self.F = bed
	88
	89
	90
[2187]	91	def tearDown(self):
	92	pass
	93
	94	def test_datapoint_at_centroid(self):
	95	a = [0.0, 0.0]
	96	b = [0.0, 2.0]
	97	c = [2.0,0.0]
	98	points = [a, b, c]
	99	vertices = [ [1,0,2] ] #bac
	100
	101	data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
	102
	103	interp = Interpolate(points, vertices)
	104	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	105	[[1./3, 1./3, 1./3]])
	106
	107	def test_quad_tree(self):
	108	p0 = [-10.0, -10.0]
	109	p1 = [20.0, -10.0]
	110	p2 = [-10.0, 20.0]
	111	p3 = [10.0, 50.0]
	112	p4 = [30.0, 30.0]
	113	p5 = [50.0, 10.0]
	114	p6 = [40.0, 60.0]
	115	p7 = [60.0, 40.0]
	116	p8 = [-66.0, 20.0]
	117	p9 = [10.0, -66.0]
	118
	119	points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
	120	triangles = [ [0, 1, 2],
	121	[3, 2, 4],
	122	[4, 2, 1],
	123	[4, 1, 5],
	124	[3, 4, 6],
	125	[6, 4, 7],
	126	[7, 4, 5],
	127	[8, 0, 2],
	128	[0, 9, 1]]
	129
	130	data = [ [4,4] ]
	131	interp = Interpolate(points, triangles,
[2201]	132	max_vertices_per_cell = 4)
[2187]	133	#print "PDSG - interp.get_A()", interp.get_A()
	134	answer = [ [ 0.06666667, 0.46666667, 0.46666667, 0.,
	135	0., 0. , 0., 0., 0., 0.]]
	136
	137
	138	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	139	answer)
	140	#interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	141	#print "PDSG - interp.get_A()", interp.get_A()
	142	data = [[-30, -30]]
	143	answer = [ [ 0.0, 0.0, 0.0, 0.,
	144	0., 0. , 0., 0., 0., 0.]]
	145
	146	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	147	answer)
	148
	149
	150	#point outside of quad tree root cell
	151	#interp.set_point_coordinates([[-70, -70]])
	152	#print "PDSG - interp.get_A()", interp.get_A()
	153	data = [[-70, -70]]
	154	answer = [ [ 0.0, 0.0, 0.0, 0.,
	155	0., 0. , 0., 0., 0., 0.]]
	156	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	157	answer)
	158
	159	def test_datapoints_at_vertices(self):
[3019]	160	#Test that data points coinciding with vertices yield a diagonal matrix
	161
[2187]	162
	163	a = [0.0, 0.0]
	164	b = [0.0, 2.0]
	165	c = [2.0,0.0]
	166	points = [a, b, c]
	167	vertices = [ [1,0,2] ] #bac
	168
	169	data = points #Use data at vertices
	170
	171	interp = Interpolate(points, vertices)
	172	answer = [[1., 0., 0.],
	173	[0., 1., 0.],
	174	[0., 0., 1.]]
	175	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	176	answer)
	177
	178
	179	def test_datapoints_on_edge_midpoints(self):
[3019]	180	#Try datapoints midway on edges -
	181	#each point should affect two matrix entries equally
	182
[2187]	183
	184	a = [0.0, 0.0]
	185	b = [0.0, 2.0]
	186	c = [2.0,0.0]
	187	points = [a, b, c]
	188	vertices = [ [1,0,2] ] #bac
	189
	190	data = [ [0., 1.], [1., 0.], [1., 1.] ]
	191	answer = [[0.5, 0.5, 0.0], #Affects vertex 1 and 0
	192	[0.5, 0.0, 0.5], #Affects vertex 0 and 2
	193	[0.0, 0.5, 0.5]]
[2690]	194	interp = Interpolate(points, vertices)
[2187]	195
	196	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	197	answer)
	198
	199	def test_datapoints_on_edges(self):
[3019]	200	#Try datapoints on edges -
	201	#each point should affect two matrix entries in proportion
	202
[2187]	203
	204	a = [0.0, 0.0]
	205	b = [0.0, 2.0]
	206	c = [2.0,0.0]
	207	points = [a, b, c]
	208	vertices = [ [1,0,2] ] #bac
	209
	210	data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
	211	answer = [[0.25, 0.75, 0.0], #Affects vertex 1 and 0
	212	[0.25, 0.0, 0.75], #Affects vertex 0 and 2
	213	[0.0, 0.25, 0.75]]
	214
[2690]	215	interp = Interpolate(points, vertices)
[2187]	216
	217	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	218	answer)
	219
	220
	221	def test_arbitrary_datapoints(self):
[3019]	222	#Try arbitrary datapoints
	223
[2187]	224
	225	from Numeric import sum
	226
	227	a = [0.0, 0.0]
	228	b = [0.0, 2.0]
	229	c = [2.0,0.0]
	230	points = [a, b, c]
	231	vertices = [ [1,0,2] ] #bac
	232
	233	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
	234
[2690]	235	interp = Interpolate(points, vertices)
[2187]	236	#print "interp.get_A()", interp.get_A()
	237	results = interp._build_interpolation_matrix_A(data).todense()
	238	assert allclose(sum(results, axis=1), 1.0)
	239
[2655]	240	def test_arbitrary_datapoints_some_outside(self):
[3019]	241	#Try arbitrary datapoints one outside the triangle.
	242	#That one should be ignored
	243
[2187]	244
	245	from Numeric import sum
	246
	247	a = [0.0, 0.0]
	248	b = [0.0, 2.0]
	249	c = [2.0,0.0]
	250	points = [a, b, c]
	251	vertices = [ [1,0,2] ] #bac
	252
	253	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
	254
[2690]	255	interp = Interpolate(points, vertices)
[2187]	256	results = interp._build_interpolation_matrix_A(data).todense()
	257	assert allclose(sum(results, axis=1), [1,1,1,0])
	258
	259
	260
	261	# this causes a memory error in scipy.sparse
	262	def test_more_triangles(self):
	263
	264	a = [-1.0, 0.0]
	265	b = [3.0, 4.0]
	266	c = [4.0,1.0]
	267	d = [-3.0, 2.0] #3
	268	e = [-1.0,-2.0]
	269	f = [1.0, -2.0] #5
	270
	271	points = [a, b, c, d,e,f]
	272	triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
	273
	274	#Data points
	275	data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
	276	interp = Interpolate(points, triangles)
	277
	278	answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d
	279	[0.5, 0.0, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	280	[0.75, 0.25, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
	281	[0.0, 0.5, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	282	[0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
	283	[1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
	284
	285
	286	A = interp._build_interpolation_matrix_A(data).todense()
	287	for i in range(A.shape[0]):
	288	for j in range(A.shape[1]):
	289	if not allclose(A[i,j], answer[i][j]):
	290	print i,j,':',A[i,j], answer[i][j]
	291
	292
[2190]	293	#results = interp._build_interpolation_matrix_A(data).todense()
[2187]	294
[2190]	295	assert allclose(A, answer)
[2690]	296
	297	def test_geo_ref(self):
	298	v0 = [0.0, 0.0]
	299	v1 = [0.0, 5.0]
	300	v2 = [5.0, 0.0]
[2187]	301
[2690]	302	vertices_absolute = [v0, v1, v2]
	303	triangles = [ [1,0,2] ] #bac
	304
	305	geo = Geo_reference(57,100, 500)
	306
	307	vertices = geo.change_points_geo_ref(vertices_absolute)
	308	#print "vertices",vertices
	309
	310	d0 = [1.0, 1.0]
	311	d1 = [1.0, 2.0]
	312	d2 = [3.0, 1.0]
	313	point_coords = [ d0, d1, d2]
	314
	315	interp = Interpolate(vertices, triangles, mesh_origin=geo)
	316	f = linear_function(vertices_absolute)
	317	z = interp.interpolate(f, point_coords)
	318	answer = linear_function(point_coords)
	319
	320	#print "z",z
	321	#print "answer",answer
	322	assert allclose(z, answer)
[2879]	323
[2690]	324
[2879]	325	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	326	answer = linear_function(point_coords)
	327
	328	#print "z",z
	329	#print "answer",answer
	330	assert allclose(z, answer)
	331
	332
[2690]	333	def test_Geospatial_verts(self):
	334	v0 = [0.0, 0.0]
	335	v1 = [0.0, 5.0]
	336	v2 = [5.0, 0.0]
	337
	338	vertices_absolute = [v0, v1, v2]
	339	triangles = [ [1,0,2] ] #bac
	340
	341	geo = Geo_reference(57,100, 500)
	342	vertices = geo.change_points_geo_ref(vertices_absolute)
	343	geopoints = Geospatial_data(vertices,geo_reference = geo)
	344	#print "vertices",vertices
	345
	346	d0 = [1.0, 1.0]
	347	d1 = [1.0, 2.0]
	348	d2 = [3.0, 1.0]
	349	point_coords = [ d0, d1, d2]
	350
	351	interp = Interpolate(geopoints, triangles)
	352	f = linear_function(vertices_absolute)
	353	z = interp.interpolate(f, point_coords)
	354	answer = linear_function(point_coords)
	355
	356	#print "z",z
	357	#print "answer",answer
	358	assert allclose(z, answer)
	359
[2879]	360	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	361	answer = linear_function(point_coords)
	362
	363	#print "z",z
	364	#print "answer",answer
	365	assert allclose(z, answer)
	366
[2189]	367	def test_interpolate_attributes_to_points(self):
	368	v0 = [0.0, 0.0]
	369	v1 = [0.0, 5.0]
	370	v2 = [5.0, 0.0]
[2187]	371
[2189]	372	vertices = [v0, v1, v2]
	373	triangles = [ [1,0,2] ] #bac
[2187]	374
[2189]	375	d0 = [1.0, 1.0]
	376	d1 = [1.0, 2.0]
	377	d2 = [3.0, 1.0]
	378	point_coords = [ d0, d1, d2]
	379
[2690]	380	interp = Interpolate(vertices, triangles)
[2189]	381	f = linear_function(vertices)
	382	z = interp.interpolate(f, point_coords)
	383	answer = linear_function(point_coords)
	384
[2690]	385	#print "z",z
	386	#print "answer",answer
[2189]	387	assert allclose(z, answer)
	388
	389
[2879]	390	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	391	answer = linear_function(point_coords)
[2189]	392
[2879]	393	#print "z",z
	394	#print "answer",answer
	395	assert allclose(z, answer)
	396
[2189]	397	def test_interpolate_attributes_to_pointsII(self):
	398	a = [-1.0, 0.0]
	399	b = [3.0, 4.0]
	400	c = [4.0, 1.0]
	401	d = [-3.0, 2.0] #3
	402	e = [-1.0, -2.0]
	403	f = [1.0, -2.0] #5
	404
	405	vertices = [a, b, c, d,e,f]
	406	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	407
	408
	409	point_coords = [[-2.0, 2.0],
	410	[-1.0, 1.0],
	411	[0.0, 2.0],
	412	[1.0, 1.0],
	413	[2.0, 1.0],
	414	[0.0, 0.0],
	415	[1.0, 0.0],
	416	[0.0, -1.0],
	417	[-0.2, -0.5],
	418	[-0.9, -1.5],
	419	[0.5, -1.9],
	420	[3.0, 1.0]]
	421
	422	interp = Interpolate(vertices, triangles)
	423	f = linear_function(vertices)
	424	z = interp.interpolate(f, point_coords)
	425	answer = linear_function(point_coords)
	426	#print "z",z
	427	#print "answer",answer
	428	assert allclose(z, answer)
	429
[2879]	430	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	431	answer = linear_function(point_coords)
	432
	433	#print "z",z
	434	#print "answer",answer
	435	assert allclose(z, answer)
	436
[2189]	437	def test_interpolate_attributes_to_pointsIII(self):
[3019]	438	#Test linear interpolation of known values at vertices to
	439	#new points inside a triangle
	440
[2189]	441	a = [0.0, 0.0]
	442	b = [0.0, 5.0]
	443	c = [5.0, 0.0]
	444	d = [5.0, 5.0]
	445
	446	vertices = [a, b, c, d]
	447	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	448
	449	#Points within triangle 1
	450	d0 = [1.0, 1.0]
	451	d1 = [1.0, 2.0]
	452	d2 = [3.0, 1.0]
	453
	454	#Point within triangle 2
	455	d3 = [4.0, 3.0]
	456
	457	#Points on common edge
	458	d4 = [2.5, 2.5]
	459	d5 = [4.0, 1.0]
	460
	461	#Point on common vertex
	462	d6 = [0., 5.]
	463
	464	point_coords = [d0, d1, d2, d3, d4, d5, d6]
	465
	466	interp = Interpolate(vertices, triangles)
	467
	468	#Known values at vertices
	469	#Functions are x+y, x+2y, 2x+y, x-y-5
	470	f = [ [0., 0., 0., -5.], # (0,0)
	471	[5., 10., 5., -10.], # (0,5)
	472	[5., 5., 10.0, 0.], # (5,0)
	473	[10., 15., 15., -5.]] # (5,5)
	474
	475	z = interp.interpolate(f, point_coords)
	476	answer = [ [2., 3., 3., -5.], # (1,1)
	477	[3., 5., 4., -6.], # (1,2)
	478	[4., 5., 7., -3.], # (3,1)
	479	[7., 10., 11., -4.], # (4,3)
	480	[5., 7.5, 7.5, -5.], # (2.5, 2.5)
	481	[5., 6., 9., -2.], # (4,1)
	482	[5., 10., 5., -10.]] # (0,5)
	483
	484	#print "***********"
	485	#print "z",z
	486	#print "answer",answer
	487	#print "***********"
	488
	489	assert allclose(z, answer)
	490
	491
[2879]	492	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	493
	494	#print "z",z
	495	#print "answer",answer
	496	assert allclose(z, answer)
	497
[2189]	498	def test_interpolate_point_outside_of_mesh(self):
[3019]	499	#Test linear interpolation of known values at vertices to
	500	#new points inside a triangle
	501
[2189]	502	a = [0.0, 0.0]
	503	b = [0.0, 5.0]
	504	c = [5.0, 0.0]
	505	d = [5.0, 5.0]
	506
	507	vertices = [a, b, c, d]
	508	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	509
	510	#Far away point
	511	d7 = [-1., -1.]
	512
	513	point_coords = [ d7]
	514	interp = Interpolate(vertices, triangles)
	515
	516	#Known values at vertices
	517	#Functions are x+y, x+2y, 2x+y, x-y-5
	518	f = [ [0., 0., 0., -5.], # (0,0)
	519	[5., 10., 5., -10.], # (0,5)
	520	[5., 5., 10.0, 0.], # (5,0)
	521	[10., 15., 15., -5.]] # (5,5)
	522
[3018]	523	z = interp.interpolate(f, point_coords) #, verbose=True)
[2655]	524	answer = array([ [INF, INF, INF, INF]]) # (-1,-1)
[2189]	525
	526	#print "***********"
	527	#print "z",z
	528	#print "answer",answer
	529	#print "***********"
	530
	531	#Should an error message be returned if points are outside
[3018]	532	# of the mesh?
	533	# A warning message is printed, if verbose is on.
[2189]	534
[2655]	535	for i in range(4):
	536	self.failUnless( z[0,i] == answer[0,i], 'Fail!')
	537
[2879]	538	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	539
	540	#print "z",z
	541	#print "answer",answer
	542
[2881]	543	for i in range(4):
	544	self.failUnless( z[0,i] == answer[0,i], 'Fail!')
	545
	546
[2189]	547	def test_interpolate_attributes_to_pointsIV(self):
	548	a = [-1.0, 0.0]
	549	b = [3.0, 4.0]
	550	c = [4.0, 1.0]
	551	d = [-3.0, 2.0] #3
	552	e = [-1.0, -2.0]
	553	f = [1.0, -2.0] #5
	554
	555	vertices = [a, b, c, d,e,f]
	556	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	557
	558
	559	point_coords = [[-2.0, 2.0],
	560	[-1.0, 1.0],
	561	[0.0, 2.0],
	562	[1.0, 1.0],
	563	[2.0, 1.0],
	564	[0.0, 0.0],
	565	[1.0, 0.0],
	566	[0.0, -1.0],
	567	[-0.2, -0.5],
	568	[-0.9, -1.5],
	569	[0.5, -1.9],
	570	[3.0, 1.0]]
	571
	572	interp = Interpolate(vertices, triangles)
	573	f = array([linear_function(vertices),2*linear_function(vertices) ])
	574	f = transpose(f)
	575	#print "f",f
	576	z = interp.interpolate(f, point_coords)
	577	answer = [linear_function(point_coords),
	578	2*linear_function(point_coords) ]
	579	answer = transpose(answer)
	580	#print "z",z
	581	#print "answer",answer
	582	assert allclose(z, answer)
	583
[2879]	584	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
[2189]	585
[2879]	586	#print "z",z
	587	#print "answer",answer
	588	assert allclose(z, answer)
	589
[2189]	590	def test_interpolate_blocking(self):
	591	a = [-1.0, 0.0]
	592	b = [3.0, 4.0]
	593	c = [4.0, 1.0]
	594	d = [-3.0, 2.0] #3
	595	e = [-1.0, -2.0]
	596	f = [1.0, -2.0] #5
	597
	598	vertices = [a, b, c, d,e,f]
	599	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	600
	601
	602	point_coords = [[-2.0, 2.0],
	603	[-1.0, 1.0],
	604	[0.0, 2.0],
	605	[1.0, 1.0],
	606	[2.0, 1.0],
	607	[0.0, 0.0],
	608	[1.0, 0.0],
	609	[0.0, -1.0],
	610	[-0.2, -0.5],
	611	[-0.9, -1.5],
	612	[0.5, -1.9],
	613	[3.0, 1.0]]
	614
	615	interp = Interpolate(vertices, triangles)
	616	f = array([linear_function(vertices),2*linear_function(vertices) ])
	617	f = transpose(f)
	618	#print "f",f
	619	for blocking_max in range(len(point_coords)+2):
	620	#if True:
	621	# blocking_max = 5
	622	z = interp.interpolate(f, point_coords,
	623	start_blocking_len=blocking_max)
	624	answer = [linear_function(point_coords),
	625	2*linear_function(point_coords) ]
	626	answer = transpose(answer)
	627	#print "z",z
	628	#print "answer",answer
	629	assert allclose(z, answer)
[2660]	630
	631	f = array([linear_function(vertices),2*linear_function(vertices),
	632	2*linear_function(vertices) - 100 ])
	633	f = transpose(f)
	634	#print "f",f
	635	for blocking_max in range(len(point_coords)+2):
	636	#if True:
	637	# blocking_max = 5
	638	z = interp.interpolate(f, point_coords,
	639	start_blocking_len=blocking_max)
	640	answer = array([linear_function(point_coords),
	641	2*linear_function(point_coords) ,
	642	2*linear_function(point_coords)-100 ])
	643	z = transpose(z)
	644	#print "z",z
	645	#print "answer",answer
	646	assert allclose(z, answer)
[2189]	647
[2787]	648	def test_interpolate_geo_spatial(self):
	649	a = [-1.0, 0.0]
	650	b = [3.0, 4.0]
	651	c = [4.0, 1.0]
	652	d = [-3.0, 2.0] #3
	653	e = [-1.0, -2.0]
	654	f = [1.0, -2.0] #5
	655
	656	vertices = [a, b, c, d,e,f]
	657	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	658
	659
	660	point_coords_absolute = [[-2.0, 2.0],
	661	[-1.0, 1.0],
	662	[0.0, 2.0],
	663	[1.0, 1.0],
	664	[2.0, 1.0],
	665	[0.0, 0.0],
	666	[1.0, 0.0],
	667	[0.0, -1.0],
	668	[-0.2, -0.5],
	669	[-0.9, -1.5],
	670	[0.5, -1.9],
	671	[3.0, 1.0]]
	672
	673	geo = Geo_reference(57,100, 500)
	674	point_coords = geo.change_points_geo_ref(point_coords_absolute)
	675	point_coords = Geospatial_data(point_coords,geo_reference = geo)
	676
	677	interp = Interpolate(vertices, triangles)
	678	f = array([linear_function(vertices),2*linear_function(vertices) ])
	679	f = transpose(f)
	680	#print "f",f
	681	for blocking_max in range(14):
	682	#if True:
	683	# blocking_max = 5
	684	z = interp.interpolate(f, point_coords,
	685	start_blocking_len=blocking_max)
	686	answer = [linear_function(point_coords.get_data_points( \
	687	absolute = True)),
	688	2*linear_function(point_coords.get_data_points( \
	689	absolute = True)) ]
	690	answer = transpose(answer)
	691	#print "z",z
	692	#print "answer",answer
	693	assert allclose(z, answer)
	694
	695	f = array([linear_function(vertices),2*linear_function(vertices),
	696	2*linear_function(vertices) - 100 ])
	697	f = transpose(f)
	698	#print "f",f
	699	for blocking_max in range(14):
	700	#if True:
	701	# blocking_max = 5
	702	z = interp.interpolate(f, point_coords,
	703	start_blocking_len=blocking_max)
	704	answer = array([linear_function(point_coords.get_data_points( \
	705	absolute = True)),
	706	2*linear_function(point_coords.get_data_points( \
	707	absolute = True)) ,
	708	2*linear_function(point_coords.get_data_points( \
	709	absolute = True))-100 ])
	710	z = transpose(z)
	711	#print "z",z
	712	#print "answer",answer
	713	assert allclose(z, answer)
	714
[2879]	715	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	716
	717	#print "z",z
	718	#print "answer",answer
	719	assert allclose(z, answer)
	720
[2787]	721	def test_interpolate_geo_spatial(self):
	722	a = [-1.0, 0.0]
	723	b = [3.0, 4.0]
	724	c = [4.0, 1.0]
	725	d = [-3.0, 2.0] #3
	726	e = [-1.0, -2.0]
	727	f = [1.0, -2.0] #5
	728
	729	vertices = [a, b, c, d,e,f]
	730	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	731
	732
	733	point_coords_absolute = [[-2.0, 2.0],
	734	[-1.0, 1.0],
	735	[0.0, 2.0],
	736	[1.0, 1.0],
	737	[2.0, 1.0],
	738	[0.0, 0.0],
	739	[1.0, 0.0],
	740	[0.0, -1.0],
	741	[-0.2, -0.5],
	742	[-0.9, -1.5],
	743	[0.5, -1.9],
	744	[3.0, 1.0]]
	745
	746	geo = Geo_reference(57,100, 500)
	747	point_coords = geo.change_points_geo_ref(point_coords_absolute)
	748	point_coords = Geospatial_data(point_coords,geo_reference = geo)
	749
	750	interp = Interpolate(vertices, triangles)
	751	f = array([linear_function(vertices),2*linear_function(vertices) ])
	752	f = transpose(f)
	753	#print "f",f
	754	z = interp.interpolate_block(f, point_coords)
	755	answer = [linear_function(point_coords.get_data_points( \
	756	absolute = True)),
	757	2*linear_function(point_coords.get_data_points( \
	758	absolute = True)) ]
	759	answer = transpose(answer)
	760	#print "z",z
	761	#print "answer",answer
	762	assert allclose(z, answer)
	763
[2879]	764	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	765
	766	#print "z",z
	767	#print "answer",answer
	768	assert allclose(z, answer)
	769
	770
[2189]	771	def test_interpolate_reuse(self):
	772	a = [-1.0, 0.0]
	773	b = [3.0, 4.0]
	774	c = [4.0, 1.0]
	775	d = [-3.0, 2.0] #3
	776	e = [-1.0, -2.0]
	777	f = [1.0, -2.0] #5
	778
	779	vertices = [a, b, c, d,e,f]
	780	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	781
	782
	783	point_coords = [[-2.0, 2.0],
	784	[-1.0, 1.0],
	785	[0.0, 2.0],
	786	[1.0, 1.0],
	787	[2.0, 1.0],
	788	[0.0, 0.0],
	789	[1.0, 0.0],
	790	[0.0, -1.0],
	791	[-0.2, -0.5],
	792	[-0.9, -1.5],
	793	[0.5, -1.9],
	794	[3.0, 1.0]]
	795
	796	interp = Interpolate(vertices, triangles)
	797	f = array([linear_function(vertices),2*linear_function(vertices) ])
	798	f = transpose(f)
	799	z = interp.interpolate(f, point_coords,
	800	start_blocking_len=20)
	801	answer = [linear_function(point_coords),
	802	2*linear_function(point_coords) ]
	803	answer = transpose(answer)
	804	#print "z",z
	805	#print "answer",answer
	806	assert allclose(z, answer)
[2201]	807	assert allclose(interp._A_can_be_reused, True)
[2189]	808
	809	z = interp.interpolate(f)
	810	assert allclose(z, answer)
	811
	812	# This causes blocking to occur.
	813	z = interp.interpolate(f, start_blocking_len=10)
	814	assert allclose(z, answer)
[2201]	815	assert allclose(interp._A_can_be_reused, False)
[2189]	816
	817	#A is recalculated
	818	z = interp.interpolate(f)
	819	assert allclose(z, answer)
[2201]	820	assert allclose(interp._A_can_be_reused, True)
	821
[2189]	822	interp = Interpolate(vertices, triangles)
	823	#Must raise an exception, no points specified
	824	try:
	825	z = interp.interpolate(f)
	826	except:
	827	pass
	828
	829
[2394]	830
	831	def test_interpolation_interface_time_only(self):
	832
[2655]	833	# Test spatio-temporal interpolation
	834	# Test that spatio temporal function performs the correct
	835	# interpolations in both time and space
	836
[2394]	837
[2655]	838
[2394]	839	#Three timesteps
	840	time = [1.0, 5.0, 6.0]
	841
	842
	843	#One quantity
	844	Q = zeros( (3,6), Float )
	845
	846	#Linear in time and space
	847	a = [0.0, 0.0]
	848	b = [0.0, 2.0]
	849	c = [2.0, 0.0]
	850	d = [0.0, 4.0]
	851	e = [2.0, 2.0]
	852	f = [4.0, 0.0]
	853
	854	points = [a, b, c, d, e, f]
	855
	856	for i, t in enumerate(time):
	857	Q[i, :] = t*linear_function(points)
	858
	859
	860	#Check basic interpolation of one quantity using averaging
	861	#(no interpolation points or spatial info)
[2563]	862	I = Interpolation_function(time, [mean(Q[0,:]),
[2394]	863	mean(Q[1,:]),
	864	mean(Q[2,:])])
	865
	866
	867
	868	#Check temporal interpolation
	869	for i in [0,1,2]:
	870	assert allclose(I(time[i]), mean(Q[i,:]))
	871
	872	#Midway
	873	assert allclose(I( (time[0] + time[1])/2 ),
	874	(I(time[0]) + I(time[1]))/2 )
	875
	876	assert allclose(I( (time[1] + time[2])/2 ),
	877	(I(time[1]) + I(time[2]))/2 )
	878
	879	assert allclose(I( (time[0] + time[2])/2 ),
	880	(I(time[0]) + I(time[2]))/2 )
	881
	882	#1/3
	883	assert allclose(I( (time[0] + time[2])/3 ),
	884	(I(time[0]) + I(time[2]))/3 )
	885
	886
	887	#Out of bounds checks
	888	try:
	889	I(time[0]-1)
	890	except:
	891	pass
	892	else:
	893	raise 'Should raise exception'
	894
	895	try:
	896	I(time[-1]+1)
	897	except:
	898	pass
	899	else:
	900	raise 'Should raise exception'
	901
	902
	903
	904
	905	def test_interpolation_interface_spatial_only(self):
[2655]	906	# Test spatio-temporal interpolation with constant time
	907
[2394]	908	#Three timesteps
	909	time = [1.0, 5.0, 6.0]
	910
	911
	912	#Setup mesh used to represent fitted function
	913	a = [0.0, 0.0]
	914	b = [0.0, 2.0]
	915	c = [2.0, 0.0]
	916	d = [0.0, 4.0]
	917	e = [2.0, 2.0]
	918	f = [4.0, 0.0]
	919
	920	points = [a, b, c, d, e, f]
	921	#bac, bce, ecf, dbe
	922	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	923
	924
	925	#New datapoints where interpolated values are sought
	926	interpolation_points = [[ 0.0, 0.0],
	927	[ 0.5, 0.5],
	928	[ 0.7, 0.7],
	929	[ 1.0, 0.5],
	930	[ 2.0, 0.4],
	931	[ 2.8, 1.2]]
	932
	933
	934	#One quantity linear in space
	935	Q = linear_function(points)
	936
	937
	938	#Check interpolation of one quantity using interpolaton points
[2563]	939	I = Interpolation_function(time, Q,
[2394]	940	vertex_coordinates = points,
	941	triangles = triangles,
	942	interpolation_points = interpolation_points,
	943	verbose = False)
	944
	945
	946	answer = linear_function(interpolation_points)
	947
	948	t = time[0]
	949	for j in range(50): #t in [1, 6]
	950	for id in range(len(interpolation_points)):
	951	assert allclose(I(t, id), answer[id])
	952
	953	t += 0.1
	954
	955
	956	try:
	957	I(1)
	958	except:
	959	pass
	960	else:
	961	raise 'Should raise exception'
	962
	963
	964
	965	def test_interpolation_interface(self):
[2655]	966	# Test spatio-temporal interpolation
	967	# Test that spatio temporal function performs the correct
	968	# interpolations in both time and space
	969
[2394]	970	#Three timesteps
[2659]	971	time = [1.0, 5.0, 6.0]
[2394]	972
	973	#Setup mesh used to represent fitted function
	974	a = [0.0, 0.0]
	975	b = [0.0, 2.0]
	976	c = [2.0, 0.0]
	977	d = [0.0, 4.0]
	978	e = [2.0, 2.0]
	979	f = [4.0, 0.0]
	980
	981	points = [a, b, c, d, e, f]
	982	#bac, bce, ecf, dbe
	983	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	984
	985
	986	#New datapoints where interpolated values are sought
	987	interpolation_points = [[ 0.0, 0.0],
	988	[ 0.5, 0.5],
	989	[ 0.7, 0.7],
	990	[ 1.0, 0.5],
	991	[ 2.0, 0.4],
	992	[ 2.8, 1.2]]
	993
	994	#One quantity
	995	Q = zeros( (3,6), Float )
	996
	997	#Linear in time and space
	998	for i, t in enumerate(time):
	999	Q[i, :] = t*linear_function(points)
	1000
	1001	#Check interpolation of one quantity using interpolaton points)
[2563]	1002	I = Interpolation_function(time, Q,
[2394]	1003	vertex_coordinates = points,
	1004	triangles = triangles,
	1005	interpolation_points = interpolation_points,
	1006	verbose = False)
	1007
	1008	answer = linear_function(interpolation_points)
	1009
	1010	t = time[0]
	1011	for j in range(50): #t in [1, 6]
	1012	for id in range(len(interpolation_points)):
	1013	assert allclose(I(t, id), t*answer[id])
	1014	t += 0.1
	1015
	1016	try:
	1017	I(1)
	1018	except:
	1019	pass
	1020	else:
	1021	raise 'Should raise exception'
	1022
	1023
[2750]	1024	def test_interpolation_precompute_points(self):
[2684]	1025	# looking at a discrete mesh
	1026	#
	1027
	1028	#Three timesteps
	1029	time = [0.0, 60.0]
	1030
	1031	#Setup mesh used to represent fitted function
	1032	points = [[ 15., -20.],
	1033	[ 15., 10.],
	1034	[ 0., -20.],
	1035	[ 0., 10.],
	1036	[ 0., -20.],
	1037	[ 15., 10.]]
	1038
	1039	triangles = [[0, 1, 2],
	1040	[3, 4, 5]]
	1041
	1042	#New datapoints where interpolated values are sought
	1043	interpolation_points = [[ 1., 0.], [0.,1.]]
	1044
	1045	#One quantity
	1046	Q = zeros( (2,6), Float )
	1047
	1048	#Linear in time and space
	1049	for i, t in enumerate(time):
	1050	Q[i, :] = t*linear_function(points)
	1051	#print "Q", Q
	1052
	1053
	1054
	1055	interp = Interpolate(points, triangles)
	1056	f = array([linear_function(points),2*linear_function(points) ])
	1057	f = transpose(f)
	1058	#print "f",f
	1059	z = interp.interpolate(f, interpolation_points)
	1060	answer = [linear_function(interpolation_points),
	1061	2*linear_function(interpolation_points) ]
	1062	answer = transpose(answer)
[2750]	1063	#print "z",z
	1064	#print "answer",answer
[2684]	1065	assert allclose(z, answer)
	1066
	1067
	1068	#Check interpolation of one quantity using interpolaton points)
	1069	I = Interpolation_function(time, Q,
	1070	vertex_coordinates = points,
	1071	triangles = triangles,
	1072	interpolation_points = interpolation_points,
	1073	verbose = False)
[2755]	1074
[2750]	1075	#print "I.precomputed_values", I.precomputed_values
[2755]	1076
	1077	msg = 'Interpolation failed'
	1078	assert allclose(I.precomputed_values['Attribute'][1], [60, 60]), msg
	1079	#self.failUnless( I.precomputed_values['Attribute'][1] == 60.0,
	1080	# ' failed')
[2684]	1081
[2659]	1082	def test_interpolation_function_outside_point(self):
	1083	# Test spatio-temporal interpolation
	1084	# Test that spatio temporal function performs the correct
	1085	# interpolations in both time and space
	1086
	1087	#Three timesteps
	1088	time = [1.0, 5.0, 6.0]
[2394]	1089
[2659]	1090	#Setup mesh used to represent fitted function
	1091	a = [0.0, 0.0]
	1092	b = [0.0, 2.0]
	1093	c = [2.0, 0.0]
	1094	d = [0.0, 4.0]
	1095	e = [2.0, 2.0]
	1096	f = [4.0, 0.0]
[2394]	1097
[2659]	1098	points = [a, b, c, d, e, f]
	1099	#bac, bce, ecf, dbe
	1100	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
[2577]	1101
	1102
[2659]	1103	#New datapoints where interpolated values are sought
	1104	interpolation_points = [[ 0.0, 0.0],
	1105	[ 0.5, 0.5],
	1106	[ 0.7, 0.7],
	1107	[ 1.0, 0.5],
	1108	[ 2.0, 0.4],
	1109	[ 545354534, 4354354353]] # outside the mesh
	1110
	1111	#One quantity
	1112	Q = zeros( (3,6), Float )
	1113
	1114	#Linear in time and space
	1115	for i, t in enumerate(time):
	1116	Q[i, :] = t*linear_function(points)
	1117
	1118	#Check interpolation of one quantity using interpolaton points)
	1119	I = Interpolation_function(time, Q,
	1120	vertex_coordinates = points,
	1121	triangles = triangles,
	1122	interpolation_points = interpolation_points,
	1123	verbose = False)
	1124
	1125	answer = linear_function(interpolation_points)
	1126
	1127	t = time[0]
	1128	for j in range(50): #t in [1, 6]
	1129	for id in range(len(interpolation_points)-1):
	1130	assert allclose(I(t, id), t*answer[id])
	1131	t += 0.1
	1132
	1133	# Now test the point outside the mesh
	1134	t = time[0]
	1135	for j in range(50): #t in [1, 6]
	1136	self.failUnless(I(t, 5) == INF, 'Fail!')
	1137	t += 0.1
	1138
	1139	try:
	1140	I(1)
	1141	except:
	1142	pass
	1143	else:
	1144	raise 'Should raise exception'
	1145
[2995]	1146
	1147	def test_interpolation_function_time(self):
[3019]	1148	#Test a long time series with an error in it (this did cause an
	1149	#error once)
	1150
[2995]	1151
	1152	time = array(\
	1153	[0.00000000e+00, 5.00000000e-02, 1.00000000e-01, 1.50000000e-01,
	1154	2.00000000e-01, 2.50000000e-01, 3.00000000e-01, 3.50000000e-01,
	1155	4.00000000e-01, 4.50000000e-01, 5.00000000e-01, 5.50000000e-01,
	1156	6.00000000e-01, 6.50000000e-01, 7.00000000e-01, 7.50000000e-01,
	1157	8.00000000e-01, 8.50000000e-01, 9.00000000e-01, 9.50000000e-01,
	1158	1.00000000e-00, 1.05000000e+00, 1.10000000e+00, 1.15000000e+00,
	1159	1.20000000e+00, 1.25000000e+00, 1.30000000e+00, 1.35000000e+00,
	1160	1.40000000e+00, 1.45000000e+00, 1.50000000e+00, 1.55000000e+00,
	1161	1.60000000e+00, 1.65000000e+00, 1.70000000e+00, 1.75000000e+00,
	1162	1.80000000e+00, 1.85000000e+00, 1.90000000e+00, 1.95000000e+00,
	1163	2.00000000e+00, 2.05000000e+00, 2.10000000e+00, 2.15000000e+00,
	1164	2.20000000e+00, 2.25000000e+00, 2.30000000e+00, 2.35000000e+00,
	1165	2.40000000e+00, 2.45000000e+00, 2.50000000e+00, 2.55000000e+00,
	1166	2.60000000e+00, 2.65000000e+00, 2.70000000e+00, 2.75000000e+00,
	1167	2.80000000e+00, 2.85000000e+00, 2.90000000e+00, 2.95000000e+00,
	1168	3.00000000e+00, 3.05000000e+00, 9.96920997e+36, 3.15000000e+00,
	1169	3.20000000e+00, 3.25000000e+00, 3.30000000e+00, 3.35000000e+00,
	1170	3.40000000e+00, 3.45000000e+00, 3.50000000e+00, 3.55000000e+00,
	1171	3.60000000e+00, 3.65000000e+00, 3.70000000e+00, 3.75000000e+00,
	1172	3.80000000e+00, 3.85000000e+00, 3.90000000e+00, 3.95000000e+00,
	1173	4.00000000e+00, 4.05000000e+00, 4.10000000e+00, 4.15000000e+00,
	1174	4.20000000e+00, 4.25000000e+00, 4.30000000e+00, 4.35000000e+00,
	1175	4.40000000e+00, 4.45000000e+00, 4.50000000e+00, 4.55000000e+00,
	1176	4.60000000e+00, 4.65000000e+00, 4.70000000e+00, 4.75000000e+00,
	1177	4.80000000e+00, 4.85000000e+00, 4.90000000e+00, 4.95000000e+00,
	1178	5.00000000e+00, 5.05000000e+00, 5.10000000e+00, 5.15000000e+00,
	1179	5.20000000e+00, 5.25000000e+00, 5.30000000e+00, 5.35000000e+00,
	1180	5.40000000e+00, 5.45000000e+00, 5.50000000e+00, 5.55000000e+00,
	1181	5.60000000e+00, 5.65000000e+00, 5.70000000e+00, 5.75000000e+00,
	1182	5.80000000e+00, 5.85000000e+00, 5.90000000e+00, 5.95000000e+00,
	1183	6.00000000e+00, 6.05000000e+00, 6.10000000e+00, 6.15000000e+00,
	1184	6.20000000e+00, 6.25000000e+00, 6.30000000e+00, 6.35000000e+00,
	1185	6.40000000e+00, 6.45000000e+00, 6.50000000e+00, 6.55000000e+00,
	1186	6.60000000e+00, 6.65000000e+00, 6.70000000e+00, 6.75000000e+00,
	1187	6.80000000e+00, 6.85000000e+00, 6.90000000e+00, 6.95000000e+00,
	1188	7.00000000e+00, 7.05000000e+00, 7.10000000e+00, 7.15000000e+00,
	1189	7.20000000e+00, 7.25000000e+00, 7.30000000e+00, 7.35000000e+00,
	1190	7.40000000e+00, 7.45000000e+00, 7.50000000e+00, 7.55000000e+00,
	1191	7.60000000e+00, 7.65000000e+00, 7.70000000e+00, 7.75000000e+00,
	1192	7.80000000e+00, 7.85000000e+00, 7.90000000e+00, 7.95000000e+00,
	1193	8.00000000e+00, 8.05000000e+00, 8.10000000e+00, 8.15000000e+00,
	1194	8.20000000e+00, 8.25000000e+00, 8.30000000e+00, 8.35000000e+00,
	1195	8.40000000e+00, 8.45000000e+00, 8.50000000e+00, 8.55000000e+00,
	1196	8.60000000e+00, 8.65000000e+00, 8.70000000e+00, 8.75000000e+00,
	1197	8.80000000e+00, 8.85000000e+00, 8.90000000e+00, 8.95000000e+00,
	1198	9.00000000e+00, 9.05000000e+00, 9.10000000e+00, 9.15000000e+00,
	1199	9.20000000e+00, 9.25000000e+00, 9.30000000e+00, 9.35000000e+00,
	1200	9.40000000e+00, 9.45000000e+00, 9.50000000e+00, 9.55000000e+00,
	1201	9.60000000e+00, 9.65000000e+00, 9.70000000e+00, 9.75000000e+00,
	1202	9.80000000e+00, 9.85000000e+00, 9.90000000e+00, 9.95000000e+00,
	1203	1.00000000e+01, 1.00500000e+01, 1.01000000e+01, 1.01500000e+01,
	1204	1.02000000e+01, 1.02500000e+01, 1.03000000e+01, 1.03500000e+01,
	1205	1.04000000e+01, 1.04500000e+01, 1.05000000e+01, 1.05500000e+01,
	1206	1.06000000e+01, 1.06500000e+01, 1.07000000e+01, 1.07500000e+01,
	1207	1.08000000e+01, 1.08500000e+01, 1.09000000e+01, 1.09500000e+01,
	1208	1.10000000e+01, 1.10500000e+01, 1.11000000e+01, 1.11500000e+01,
	1209	1.12000000e+01, 1.12500000e+01, 1.13000000e+01, 1.13500000e+01,
	1210	1.14000000e+01, 1.14500000e+01, 1.15000000e+01, 1.15500000e+01,
	1211	1.16000000e+01, 1.16500000e+01, 1.17000000e+01, 1.17500000e+01,
	1212	1.18000000e+01, 1.18500000e+01, 1.19000000e+01, 1.19500000e+01,
	1213	1.20000000e+01, 1.20500000e+01, 1.21000000e+01, 1.21500000e+01,
	1214	1.22000000e+01, 1.22500000e+01, 1.23000000e+01, 1.23500000e+01,
	1215	1.24000000e+01, 1.24500000e+01, 1.25000000e+01, 1.25500000e+01,
	1216	1.26000000e+01, 1.26500000e+01, 1.27000000e+01, 1.27500000e+01,
	1217	1.28000000e+01, 1.28500000e+01, 1.29000000e+01, 1.29500000e+01,
	1218	1.30000000e+01, 1.30500000e+01, 1.31000000e+01, 1.31500000e+01,
	1219	1.32000000e+01, 1.32500000e+01, 1.33000000e+01, 1.33500000e+01,
	1220	1.34000000e+01, 1.34500000e+01, 1.35000000e+01, 1.35500000e+01,
	1221	1.36000000e+01, 1.36500000e+01, 1.37000000e+01, 1.37500000e+01,
	1222	1.38000000e+01, 1.38500000e+01, 1.39000000e+01, 1.39500000e+01,
	1223	1.40000000e+01, 1.40500000e+01, 1.41000000e+01, 1.41500000e+01,
	1224	1.42000000e+01, 1.42500000e+01, 1.43000000e+01, 1.43500000e+01,
	1225	1.44000000e+01, 1.44500000e+01, 1.45000000e+01, 1.45500000e+01,
	1226	1.46000000e+01, 1.46500000e+01, 1.47000000e+01, 1.47500000e+01,
	1227	1.48000000e+01, 1.48500000e+01, 1.49000000e+01, 1.49500000e+01,
	1228	1.50000000e+01, 1.50500000e+01, 1.51000000e+01, 1.51500000e+01,
	1229	1.52000000e+01, 1.52500000e+01, 1.53000000e+01, 1.53500000e+01,
	1230	1.54000000e+01, 1.54500000e+01, 1.55000000e+01, 1.55500000e+01,
	1231	1.56000000e+01, 1.56500000e+01, 1.57000000e+01, 1.57500000e+01,
	1232	1.58000000e+01, 1.58500000e+01, 1.59000000e+01, 1.59500000e+01,
	1233	1.60000000e+01, 1.60500000e+01, 1.61000000e+01, 1.61500000e+01,
	1234	1.62000000e+01, 1.62500000e+01, 1.63000000e+01, 1.63500000e+01,
	1235	1.64000000e+01, 1.64500000e+01, 1.65000000e+01, 1.65500000e+01,
	1236	1.66000000e+01, 1.66500000e+01, 1.67000000e+01, 1.67500000e+01,
	1237	1.68000000e+01, 1.68500000e+01, 1.69000000e+01, 1.69500000e+01,
	1238	1.70000000e+01, 1.70500000e+01, 1.71000000e+01, 1.71500000e+01,
	1239	1.72000000e+01, 1.72500000e+01, 1.73000000e+01, 1.73500000e+01,
	1240	1.74000000e+01, 1.74500000e+01, 1.75000000e+01, 1.75500000e+01,
	1241	1.76000000e+01, 1.76500000e+01, 1.77000000e+01, 1.77500000e+01,
	1242	1.78000000e+01, 1.78500000e+01, 1.79000000e+01, 1.79500000e+01,
	1243	1.80000000e+01, 1.80500000e+01, 1.81000000e+01, 1.81500000e+01,
	1244	1.82000000e+01, 1.82500000e+01, 1.83000000e+01, 1.83500000e+01,
	1245	1.84000000e+01, 1.84500000e+01, 1.85000000e+01, 1.85500000e+01,
	1246	1.86000000e+01, 1.86500000e+01, 1.87000000e+01, 1.87500000e+01,
	1247	1.88000000e+01, 1.88500000e+01, 1.89000000e+01, 1.89500000e+01,
	1248	1.90000000e+01, 1.90500000e+01, 1.91000000e+01, 1.91500000e+01,
	1249	1.92000000e+01, 1.92500000e+01, 1.93000000e+01, 1.93500000e+01,
	1250	1.94000000e+01, 1.94500000e+01, 1.95000000e+01, 1.95500000e+01,
	1251	1.96000000e+01, 1.96500000e+01, 1.97000000e+01, 1.97500000e+01,
	1252	1.98000000e+01, 1.98500000e+01, 1.99000000e+01, 1.99500000e+01,
	1253	2.00000000e+01, 2.00500000e+01, 2.01000000e+01, 2.01500000e+01,
	1254	2.02000000e+01, 2.02500000e+01, 2.03000000e+01, 2.03500000e+01,
	1255	2.04000000e+01, 2.04500000e+01, 2.05000000e+01, 2.05500000e+01,
	1256	2.06000000e+01, 2.06500000e+01, 2.07000000e+01, 2.07500000e+01,
	1257	2.08000000e+01, 2.08500000e+01, 2.09000000e+01, 2.09500000e+01,
	1258	2.10000000e+01, 2.10500000e+01, 2.11000000e+01, 2.11500000e+01,
	1259	2.12000000e+01, 2.12500000e+01, 2.13000000e+01, 2.13500000e+01,
	1260	2.14000000e+01, 2.14500000e+01, 2.15000000e+01, 2.15500000e+01,
	1261	2.16000000e+01, 2.16500000e+01, 2.17000000e+01, 2.17500000e+01,
	1262	2.18000000e+01, 2.18500000e+01, 2.19000000e+01, 2.19500000e+01,
	1263	2.20000000e+01, 2.20500000e+01, 2.21000000e+01, 2.21500000e+01,
	1264	2.22000000e+01, 2.22500000e+01, 2.23000000e+01, 2.23500000e+01,
	1265	2.24000000e+01, 2.24500000e+01, 2.25000000e+01])
	1266
	1267	#print 'Diff', time[1:] - time[:-1]
	1268
	1269	#Setup mesh used to represent fitted function
	1270	a = [0.0, 0.0]
	1271	b = [0.0, 2.0]
	1272	c = [2.0, 0.0]
	1273	d = [0.0, 4.0]
	1274	e = [2.0, 2.0]
	1275	f = [4.0, 0.0]
	1276
	1277	points = [a, b, c, d, e, f]
	1278	#bac, bce, ecf, dbe
	1279	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1280
	1281
	1282	#New datapoints where interpolated values are sought
	1283	interpolation_points = [[ 0.0, 0.0],
	1284	[ 0.5, 0.5],
	1285	[ 0.7, 0.7],
	1286	[ 1.0, 0.5],
	1287	[ 2.0, 0.4],
	1288	[ 545354534, 4354354353]] # outside the mesh
	1289
	1290	#One quantity
	1291	Q = zeros( (len(time),6), Float )
	1292
	1293	#Linear in time and space
	1294	for i, t in enumerate(time):
	1295	Q[i, :] = t*linear_function(points)
	1296
	1297	#Check interpolation of one quantity using interpolaton points)
	1298	try:
	1299	I = Interpolation_function(time, Q,
	1300	vertex_coordinates = points,
	1301	triangles = triangles,
	1302	interpolation_points = interpolation_points,
	1303	verbose = False)
	1304	except:
	1305	pass
	1306	else:
	1307	raise 'Should raise exception due to time being non-monotoneous'
	1308
	1309
[2655]	1310	def test_points_outside_the_polygon(self):
	1311	a = [-1.0, 0.0]
	1312	b = [3.0, 4.0]
	1313	c = [4.0, 1.0]
	1314	d = [-3.0, 2.0] #3
	1315	e = [-1.0, -2.0]
	1316	f = [1.0, -2.0] #5
	1317
	1318	vertices = [a, b, c, d,e,f]
	1319	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	1320
	1321	point_coords = [[-2.0, 2.0],
	1322	[-1.0, 1.0],
	1323	[9999.0, 9999.0], # point Outside poly
	1324	[-9999.0, 1.0], # point Outside poly
	1325	[2.0, 1.0],
	1326	[0.0, 0.0],
	1327	[1.0, 0.0],
	1328	[0.0, -1.0],
	1329	[-0.2, -0.5],
	1330	[-0.9, -1.5],
	1331	[0.5, -1.9],
	1332	[999999, 9999999]] # point Outside poly
	1333	geo_data = Geospatial_data(data_points = point_coords)
	1334
	1335	interp = Interpolate(vertices, triangles)
	1336	f = array([linear_function(vertices),2*linear_function(vertices) ])
	1337	f = transpose(f)
	1338	#print "f",f
	1339	z = interp.interpolate(f, geo_data)
	1340	#z = interp.interpolate(f, point_coords)
	1341	answer = [linear_function(point_coords),
	1342	2*linear_function(point_coords) ]
	1343	answer = transpose(answer)
	1344	answer[2,:] = [INF, INF]
	1345	answer[3,:] = [INF, INF]
	1346	answer[11,:] = [INF, INF]
	1347	#print "z",z
	1348	#print "answer _ fixed",answer
	1349	assert allclose(z[0:1], answer[0:1])
	1350	assert allclose(z[4:10], answer[4:10])
	1351	for i in [2,3,11]:
	1352	self.failUnless( z[i,1] == answer[11,1], 'Fail!')
	1353	self.failUnless( z[i,0] == answer[11,0], 'Fail!')
[3411]	1354
	1355	def ztest_interpolate_sww2csv(self):
	1356
	1357	def elevation_function(x, y):
	1358	return -x
[2655]	1359
[3411]	1360	# create mesh
	1361	mesh_file = tempfile.mktemp(".tsh")
	1362	points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
	1363	m = Mesh()
	1364	m.add_vertices(points)
	1365	m.autoSegment()
	1366	m.generate_mesh(verbose=False)
	1367	m.export_mesh_file(mesh_file)
	1368
	1369	#Create shallow water domain
	1370	domain = Domain(mesh_file)
	1371	os.remove(mesh_file)
	1372
	1373	domain.default_order=2
	1374	domain.beta_h = 0
	1375
	1376	#Set some field values
	1377	domain.set_quantity('elevation', elevation_function)
	1378	domain.set_quantity('friction', 0.03)
	1379	domain.set_quantity('xmomentum', 10.0)
	1380	domain.set_quantity('ymomentum', 100.0)
	1381
	1382	######################
	1383	# Boundary conditions
	1384	B = Transmissive_boundary(domain)
	1385	domain.set_boundary( {'exterior': B})
	1386
	1387	# This call mangles the stage values.
	1388	domain.distribute_to_vertices_and_edges()
	1389	domain.set_quantity('stage', 1.0)
	1390
	1391	#sww_file = tempfile.mktemp("")
	1392	domain.filename = 'datatest' + str(time.time())
	1393	#domain.filename = sww_file
	1394	#print "domain.filename",domain.filename
	1395	domain.format = 'sww'
	1396	domain.smooth = True
	1397	domain.reduction = mean
	1398
	1399	sww = get_dataobject(domain)
	1400	sww.store_connectivity()
	1401	sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
	1402	domain.set_quantity('stage', -1.0)
	1403	domain.time = 2.
	1404	sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
	1405
	1406	# test the function
	1407	points = [[5.0,1.],[0.5,2.]]
	1408	interpolate_sww2csv(sww.filename, points, "aa.csv", "bb.csv")
	1409
	1410	print "sww.filename",sww.filename
	1411	#os.remove(sww.filename)
[2187]	1412	#-------------------------------------------------------------
	1413	if __name__ == "__main__":
[2655]	1414
[3412]	1415	suite = unittest.makeSuite(Test_Interpolate,'test')
	1416	#suite = unittest.makeSuite(Test_Interpolate,'test_interpolate_sww2csv')
[2187]	1417	runner = unittest.TextTestRunner(verbosity=1)
	1418	runner.run(suite)
	1419
	1420
	1421
	1422
	1423

Note: See TracBrowser for help on using the repository browser.

Download in other formats: