Context Navigation

test_interpolate.py @ 6982

Last change on this file since 6982 was 6689, checked in by rwilson, 16 years ago
Back-merge from Numeric trunk to numpy branch.
File size: 64.5 KB

Rev	Line
[6174]	1	#!/usr/bin/env python
	2
	3	#TEST
	4
	5	#import time, os
	6
	7
	8	import sys
	9	import os
	10	import unittest
	11	from math import sqrt
	12	import tempfile
	13	import csv
	14
	15	from Scientific.IO.NetCDF import NetCDFFile
	16
[6304]	17	import numpy as num
[6174]	18
	19
	20
	21	# ANUGA code imports
	22	from interpolate import *
	23	from anuga.coordinate_transforms.geo_reference import Geo_reference
	24	from anuga.shallow_water import Domain, Transmissive_boundary
	25	from anuga.utilities.numerical_tools import mean, NAN
	26	from anuga.shallow_water.data_manager import get_dataobject
	27	from anuga.geospatial_data.geospatial_data import Geospatial_data
	28	from anuga.pmesh.mesh import Mesh
	29
	30	def distance(x, y):
	31	return sqrt(num.sum((num.array(x)-num.array(y))**2))
	32
	33	def linear_function(point):
	34	point = num.array(point)
	35	return point[:,0]+point[:,1]
	36
	37
	38	class Test_Interpolate(unittest.TestCase):
	39
	40	def setUp(self):
	41
	42	import time
	43	from mesh_factory import rectangular
	44
	45
	46	#Create basic mesh
	47	points, vertices, boundary = rectangular(2, 2)
	48
	49	#Create shallow water domain
	50	domain = Domain(points, vertices, boundary)
	51	domain.default_order=2
	52
	53
	54	#Set some field values
	55	domain.set_quantity('elevation', lambda x,y: -x)
	56	domain.set_quantity('friction', 0.03)
	57
	58
	59	######################
	60	# Boundary conditions
	61	B = Transmissive_boundary(domain)
	62	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
	63
	64
	65	######################
	66	#Initial condition - with jumps
	67
	68	bed = domain.quantities['elevation'].vertex_values
[6304]	69	stage = num.zeros(bed.shape, num.float)
[6174]	70
	71	h = 0.3
	72	for i in range(stage.shape[0]):
	73	if i % 2 == 0:
	74	stage[i,:] = bed[i,:] + h
	75	else:
	76	stage[i,:] = bed[i,:]
	77
	78	domain.set_quantity('stage', stage)
	79
	80	domain.distribute_to_vertices_and_edges()
	81
	82
	83	self.domain = domain
	84
	85	C = domain.get_vertex_coordinates()
	86	self.X = C[:,0:6:2].copy()
	87	self.Y = C[:,1:6:2].copy()
	88
	89	self.F = bed
	90
	91
	92
	93	def tearDown(self):
	94	pass
	95
	96	def test_datapoint_at_centroid(self):
	97	a = [0.0, 0.0]
	98	b = [0.0, 2.0]
	99	c = [2.0,0.0]
	100	points = [a, b, c]
	101	vertices = [ [1,0,2] ] #bac
	102
	103	data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
	104
	105	interp = Interpolate(points, vertices)
	106	A, _, _ = interp._build_interpolation_matrix_A(data)
	107	assert num.allclose(A.todense(), [[1./3, 1./3, 1./3]])
	108
	109
	110
	111	def test_simple_interpolation_example(self):
	112
	113	from mesh_factory import rectangular
	114	from shallow_water import Domain
	115	from abstract_2d_finite_volumes.quantity import Quantity
	116
	117	# Create basic mesh
	118	points, vertices, boundary = rectangular(1, 3)
	119
	120	# Create shallow water domain
	121	domain = Domain(points, vertices, boundary)
	122
	123	#----------------
	124	#Constant values
	125	#----------------
	126	quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
	127	[4,4,4],[5,5,5]])
	128
	129
	130	x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
[6304]	131	vertex_coordinates = num.concatenate( (x[:, num.newaxis], y[:, num.newaxis]), axis=1 )
[6174]	132	# FIXME: This concat should roll into get_vertex_values
	133
	134
	135	# Get interpolated values at centroids
	136	interpolation_points = domain.get_centroid_coordinates()
	137	answer = quantity.get_values(location='centroids')
	138
	139	I = Interpolate(vertex_coordinates, triangles)
	140	result = I.interpolate(vertex_values, interpolation_points)
	141	assert num.allclose(result, answer)
	142
	143
	144	#----------------
	145	# Variable values
	146	#----------------
	147	quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
	148	[1,4,-9],[2,5,0]])
	149
	150	x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
[6304]	151	vertex_coordinates = num.concatenate( (x[:, num.newaxis], y[:, num.newaxis]), axis=1 )
[6174]	152	# FIXME: This concat should roll into get_vertex_values
	153
	154
	155	# Get interpolated values at centroids
	156	interpolation_points = domain.get_centroid_coordinates()
	157	answer = quantity.get_values(location='centroids')
	158
	159	I = Interpolate(vertex_coordinates, triangles)
	160	result = I.interpolate(vertex_values, interpolation_points)
	161	assert num.allclose(result, answer)
	162
	163
	164	def test_simple_interpolation_example_using_direct_interface(self):
	165
	166	from mesh_factory import rectangular
	167	from shallow_water import Domain
	168	from abstract_2d_finite_volumes.quantity import Quantity
	169
	170	# Create basic mesh
	171	points, vertices, boundary = rectangular(1, 3)
	172
	173	# Create shallow water domain
	174	domain = Domain(points, vertices, boundary)
	175
	176	#----------------
	177	# Constant values
	178	#----------------
	179	quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
	180	[4,4,4],[5,5,5]])
	181
	182
	183	x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
[6304]	184	vertex_coordinates = num.concatenate( (x[:, num.newaxis], y[:, num.newaxis]), axis=1 )
[6174]	185	# FIXME: This concat should roll into get_vertex_values
	186
	187
	188	# Get interpolated values at centroids
	189	interpolation_points = domain.get_centroid_coordinates()
	190	answer = quantity.get_values(location='centroids')
	191
	192	result = interpolate(vertex_coordinates, triangles, vertex_values, interpolation_points)
	193	assert num.allclose(result, answer)
	194
	195
	196	#----------------
	197	# Variable values
	198	#----------------
	199	quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
	200	[1,4,-9],[2,5,0]])
	201
	202	x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
[6304]	203	vertex_coordinates = num.concatenate( (x[:, num.newaxis], y[:, num.newaxis]), axis=1 )
[6174]	204	# FIXME: This concat should roll into get_vertex_values
	205
	206
	207	# Get interpolated values at centroids
	208	interpolation_points = domain.get_centroid_coordinates()
	209	answer = quantity.get_values(location='centroids')
	210
	211	result = interpolate(vertex_coordinates, triangles,
	212	vertex_values, interpolation_points)
	213	assert num.allclose(result, answer)
	214
	215
	216	def test_simple_interpolation_example_using_direct_interface_and_caching(self):
	217
	218	from mesh_factory import rectangular
	219	from shallow_water import Domain
	220	from abstract_2d_finite_volumes.quantity import Quantity
	221
	222	# Create basic mesh
	223	points, vertices, boundary = rectangular(1, 3)
	224
	225	# Create shallow water domain
	226	domain = Domain(points, vertices, boundary)
	227
	228	#----------------
	229	# First call
	230	#----------------
	231	quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
	232	[1,4,-9],[2,5,0]])
	233
	234	x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
[6304]	235	vertex_coordinates = num.concatenate( (x[:, num.newaxis], y[:, num.newaxis]), axis=1 )
[6174]	236	# FIXME: This concat should roll into get_vertex_values
	237
	238
	239	# Get interpolated values at centroids
	240	interpolation_points = domain.get_centroid_coordinates()
	241	answer = quantity.get_values(location='centroids')
	242
	243	result = interpolate(vertex_coordinates, triangles,
	244	vertex_values, interpolation_points,
	245	use_cache=True,
	246	verbose=False)
	247	assert num.allclose(result, answer)
	248
	249	# Second call using the cache
	250	result = interpolate(vertex_coordinates, triangles,
	251	vertex_values, interpolation_points,
	252	use_cache=True,
	253	verbose=False)
	254	assert num.allclose(result, answer)
	255
	256
	257	def test_quad_tree(self):
	258	p0 = [-10.0, -10.0]
	259	p1 = [20.0, -10.0]
	260	p2 = [-10.0, 20.0]
	261	p3 = [10.0, 50.0]
	262	p4 = [30.0, 30.0]
	263	p5 = [50.0, 10.0]
	264	p6 = [40.0, 60.0]
	265	p7 = [60.0, 40.0]
	266	p8 = [-66.0, 20.0]
	267	p9 = [10.0, -66.0]
	268
	269	points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
	270	triangles = [ [0, 1, 2],
	271	[3, 2, 4],
	272	[4, 2, 1],
	273	[4, 1, 5],
	274	[3, 4, 6],
	275	[6, 4, 7],
	276	[7, 4, 5],
	277	[8, 0, 2],
	278	[0, 9, 1]]
	279
	280	data = [ [4,4] ]
	281	interp = Interpolate(points, triangles,
	282	max_vertices_per_cell = 4)
	283	#print "PDSG - interp.get_A()", interp.get_A()
	284	answer = [ [ 0.06666667, 0.46666667, 0.46666667, 0.,
	285	0., 0. , 0., 0., 0., 0.]]
	286
	287	A,_,_ = interp._build_interpolation_matrix_A(data)
	288	assert num.allclose(A.todense(), answer)
	289
	290	#interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	291	#print "PDSG - interp.get_A()", interp.get_A()
	292	data = [[-30, -30]]
	293	answer = [ [ 0.0, 0.0, 0.0, 0.,
	294	0., 0. , 0., 0., 0., 0.]]
	295
	296	A,_,_ = interp._build_interpolation_matrix_A(data)
	297	assert num.allclose(A.todense(), answer)
	298
	299
	300	#point outside of quad tree root cell
	301	#interp.set_point_coordinates([[-70, -70]])
	302	#print "PDSG - interp.get_A()", interp.get_A()
	303	data = [[-70, -70]]
	304	answer = [ [ 0.0, 0.0, 0.0, 0.,
	305	0., 0. , 0., 0., 0., 0.]]
	306
	307	A,_,_ = interp._build_interpolation_matrix_A(data)
	308	assert num.allclose(A.todense(), answer)
	309
	310
	311	def test_datapoints_at_vertices(self):
	312	#Test that data points coinciding with vertices yield a diagonal matrix
	313
	314
	315	a = [0.0, 0.0]
	316	b = [0.0, 2.0]
	317	c = [2.0,0.0]
	318	points = [a, b, c]
	319	vertices = [ [1,0,2] ] #bac
	320
	321	data = points #Use data at vertices
	322
	323	interp = Interpolate(points, vertices)
	324	answer = [[1., 0., 0.],
	325	[0., 1., 0.],
	326	[0., 0., 1.]]
	327
	328	A,_,_ = interp._build_interpolation_matrix_A(data)
	329	assert num.allclose(A.todense(), answer)
	330
	331
	332	def test_datapoints_on_edge_midpoints(self):
	333	#Try datapoints midway on edges -
	334	#each point should affect two matrix entries equally
	335
	336
	337	a = [0.0, 0.0]
	338	b = [0.0, 2.0]
	339	c = [2.0,0.0]
	340	points = [a, b, c]
	341	vertices = [ [1,0,2] ] #bac
	342
	343	data = [ [0., 1.], [1., 0.], [1., 1.] ]
	344	answer = [[0.5, 0.5, 0.0], #Affects vertex 1 and 0
	345	[0.5, 0.0, 0.5], #Affects vertex 0 and 2
	346	[0.0, 0.5, 0.5]]
	347	interp = Interpolate(points, vertices)
	348
	349	A,_,_ = interp._build_interpolation_matrix_A(data)
	350	assert num.allclose(A.todense(), answer)
	351
	352	def test_datapoints_on_edges(self):
	353	#Try datapoints on edges -
	354	#each point should affect two matrix entries in proportion
	355
	356
	357	a = [0.0, 0.0]
	358	b = [0.0, 2.0]
	359	c = [2.0,0.0]
	360	points = [a, b, c]
	361	vertices = [ [1,0,2] ] #bac
	362
	363	data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
	364	answer = [[0.25, 0.75, 0.0], #Affects vertex 1 and 0
	365	[0.25, 0.0, 0.75], #Affects vertex 0 and 2
	366	[0.0, 0.25, 0.75]]
	367
	368	interp = Interpolate(points, vertices)
	369
	370	A,_,_ = interp._build_interpolation_matrix_A(data)
	371	assert num.allclose(A.todense(), answer)
	372
	373
	374	def test_arbitrary_datapoints(self):
	375	#Try arbitrary datapoints
	376
	377
	378	a = [0.0, 0.0]
	379	b = [0.0, 2.0]
	380	c = [2.0,0.0]
	381	points = [a, b, c]
	382	vertices = [ [1,0,2] ] #bac
	383
	384	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
	385
	386	interp = Interpolate(points, vertices)
	387	#print "interp.get_A()", interp.get_A()
	388
	389	A,_,_ = interp._build_interpolation_matrix_A(data)
	390	results = A.todense()
	391	assert num.allclose(num.sum(results, axis=1), 1.0)
	392
	393	def test_arbitrary_datapoints_some_outside(self):
	394	#Try arbitrary datapoints one outside the triangle.
	395	#That one should be ignored
	396
	397
	398	a = [0.0, 0.0]
	399	b = [0.0, 2.0]
	400	c = [2.0,0.0]
	401	points = [a, b, c]
	402	vertices = [ [1,0,2] ] #bac
	403
	404	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
	405
	406	interp = Interpolate(points, vertices)
	407
	408	A,_,_ = interp._build_interpolation_matrix_A(data)
	409	results = A.todense()
	410	assert num.allclose(num.sum(results, axis=1), [1,1,1,0])
	411
	412
	413
	414	# this causes a memory error in scipy.sparse
	415	def test_more_triangles(self):
	416
	417	a = [-1.0, 0.0]
	418	b = [3.0, 4.0]
	419	c = [4.0,1.0]
	420	d = [-3.0, 2.0] #3
	421	e = [-1.0,-2.0]
	422	f = [1.0, -2.0] #5
	423
	424	points = [a, b, c, d,e,f]
	425	triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
	426
	427	#Data points
	428	data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
	429	interp = Interpolate(points, triangles)
	430
	431	answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d
	432	[0.5, 0.0, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	433	[0.75, 0.25, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
	434	[0.0, 0.5, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	435	[0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
	436	[1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
	437
	438
	439	A,_,_ = interp._build_interpolation_matrix_A(data)
	440	A = A.todense()
	441	for i in range(A.shape[0]):
	442	for j in range(A.shape[1]):
	443	if not num.allclose(A[i,j], answer[i][j]):
	444	print i,j,':',A[i,j], answer[i][j]
	445
	446
	447	#results = interp._build_interpolation_matrix_A(data).todense()
	448
	449	assert num.allclose(A, answer)
	450
	451	def test_geo_ref(self):
	452	v0 = [0.0, 0.0]
	453	v1 = [0.0, 5.0]
	454	v2 = [5.0, 0.0]
	455
	456	vertices_absolute = [v0, v1, v2]
	457	triangles = [ [1,0,2] ] #bac
	458
	459	geo = Geo_reference(57,100, 500)
	460
	461	vertices = geo.change_points_geo_ref(vertices_absolute)
	462	#print "vertices",vertices
	463
	464	d0 = [1.0, 1.0]
	465	d1 = [1.0, 2.0]
	466	d2 = [3.0, 1.0]
	467	point_coords = [ d0, d1, d2]
	468
	469	interp = Interpolate(vertices, triangles, mesh_origin=geo)
	470	f = linear_function(vertices_absolute)
	471	z = interp.interpolate(f, point_coords)
	472	answer = linear_function(point_coords)
	473
	474	#print "z",z
	475	#print "answer",answer
	476	assert num.allclose(z, answer)
	477
	478
	479	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	480	answer = linear_function(point_coords)
	481
	482	#print "z",z
	483	#print "answer",answer
	484	assert num.allclose(z, answer)
	485
	486
	487	def test_sigma_epsilon(self):
	488	"""
	489	def test_sigma_epsilon(self):
	490	Testing ticket 168. I could not reduce the bug to this small
	491	test though.
	492
	493	"""
	494	v0 = [22031.25, 59687.5]
	495	v1 = [22500., 60000.]
	496	v2 = [22350.31640625, 59716.71484375]
	497
	498	vertices = [v0, v1, v2]
	499	triangles = [ [1,0,2] ] #bac
	500
	501
	502	point_coords = [[22050., 59700.]]
	503
	504	interp = Interpolate(vertices, triangles)
	505	f = linear_function(vertices)
	506	z = interp.interpolate(f, point_coords)
	507	answer = linear_function(point_coords)
	508
	509	#print "z",z
	510	#print "answer",answer
	511	assert num.allclose(z, answer)
	512
	513
	514	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	515	answer = linear_function(point_coords)
	516
	517	#print "z",z
	518	#print "answer",answer
	519	assert num.allclose(z, answer)
	520
	521
	522	def test_Geospatial_verts(self):
	523	v0 = [0.0, 0.0]
	524	v1 = [0.0, 5.0]
	525	v2 = [5.0, 0.0]
	526
	527	vertices_absolute = [v0, v1, v2]
	528	triangles = [ [1,0,2] ] #bac
	529
	530	geo = Geo_reference(57,100, 500)
	531	vertices = geo.change_points_geo_ref(vertices_absolute)
	532	geopoints = Geospatial_data(vertices,geo_reference = geo)
	533	#print "vertices",vertices
	534
	535	d0 = [1.0, 1.0]
	536	d1 = [1.0, 2.0]
	537	d2 = [3.0, 1.0]
	538	point_coords = [ d0, d1, d2]
	539
	540	interp = Interpolate(geopoints, triangles)
	541	f = linear_function(vertices_absolute)
	542	z = interp.interpolate(f, point_coords)
	543	answer = linear_function(point_coords)
	544
	545	#print "z",z
	546	#print "answer",answer
	547	assert num.allclose(z, answer)
	548
	549	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	550	answer = linear_function(point_coords)
	551
	552	#print "z",z
	553	#print "answer",answer
	554	assert num.allclose(z, answer)
	555
	556	def test_interpolate_attributes_to_points(self):
	557	v0 = [0.0, 0.0]
	558	v1 = [0.0, 5.0]
	559	v2 = [5.0, 0.0]
	560
	561	vertices = [v0, v1, v2]
	562	triangles = [ [1,0,2] ] #bac
	563
	564	d0 = [1.0, 1.0]
	565	d1 = [1.0, 2.0]
	566	d2 = [3.0, 1.0]
	567	point_coords = [ d0, d1, d2]
	568
	569	interp = Interpolate(vertices, triangles)
	570	f = linear_function(vertices)
	571	z = interp.interpolate(f, point_coords)
	572	answer = linear_function(point_coords)
	573
	574	#print "z",z
	575	#print "answer",answer
	576	assert num.allclose(z, answer)
	577
	578
	579	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	580	answer = linear_function(point_coords)
	581
	582	#print "z",z
	583	#print "answer",answer
	584	assert num.allclose(z, answer)
	585
	586	def test_interpolate_attributes_to_pointsII(self):
	587	a = [-1.0, 0.0]
	588	b = [3.0, 4.0]
	589	c = [4.0, 1.0]
	590	d = [-3.0, 2.0] #3
	591	e = [-1.0, -2.0]
	592	f = [1.0, -2.0] #5
	593
	594	vertices = [a, b, c, d,e,f]
	595	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	596
	597
	598	point_coords = [[-2.0, 2.0],
	599	[-1.0, 1.0],
	600	[0.0, 2.0],
	601	[1.0, 1.0],
	602	[2.0, 1.0],
	603	[0.0, 0.0],
	604	[1.0, 0.0],
	605	[0.0, -1.0],
	606	[-0.2, -0.5],
	607	[-0.9, -1.5],
	608	[0.5, -1.9],
	609	[3.0, 1.0]]
	610
	611	interp = Interpolate(vertices, triangles)
	612	f = linear_function(vertices)
	613	z = interp.interpolate(f, point_coords)
	614	answer = linear_function(point_coords)
	615	#print "z",z
	616	#print "answer",answer
	617	assert num.allclose(z, answer)
	618
	619	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	620	answer = linear_function(point_coords)
	621
	622	#print "z",z
	623	#print "answer",answer
	624	assert num.allclose(z, answer)
	625
	626	def test_interpolate_attributes_to_pointsIII(self):
	627	#Test linear interpolation of known values at vertices to
	628	#new points inside a triangle
	629
	630	a = [0.0, 0.0]
	631	b = [0.0, 5.0]
	632	c = [5.0, 0.0]
	633	d = [5.0, 5.0]
	634
	635	vertices = [a, b, c, d]
	636	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	637
	638	#Points within triangle 1
	639	d0 = [1.0, 1.0]
	640	d1 = [1.0, 2.0]
	641	d2 = [3.0, 1.0]
	642
	643	#Point within triangle 2
	644	d3 = [4.0, 3.0]
	645
	646	#Points on common edge
	647	d4 = [2.5, 2.5]
	648	d5 = [4.0, 1.0]
	649
	650	#Point on common vertex
	651	d6 = [0., 5.]
	652
	653	point_coords = [d0, d1, d2, d3, d4, d5, d6]
	654
	655	interp = Interpolate(vertices, triangles)
	656
	657	#Known values at vertices
	658	#Functions are x+y, x+2y, 2x+y, x-y-5
	659	f = [ [0., 0., 0., -5.], # (0,0)
	660	[5., 10., 5., -10.], # (0,5)
	661	[5., 5., 10.0, 0.], # (5,0)
	662	[10., 15., 15., -5.]] # (5,5)
	663
	664	z = interp.interpolate(f, point_coords)
	665	answer = [ [2., 3., 3., -5.], # (1,1)
	666	[3., 5., 4., -6.], # (1,2)
	667	[4., 5., 7., -3.], # (3,1)
	668	[7., 10., 11., -4.], # (4,3)
	669	[5., 7.5, 7.5, -5.], # (2.5, 2.5)
	670	[5., 6., 9., -2.], # (4,1)
	671	[5., 10., 5., -10.]] # (0,5)
	672
	673	#print "***********"
	674	#print "z",z
	675	#print "answer",answer
	676	#print "***********"
	677
	678	assert num.allclose(z, answer)
	679
	680
	681	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	682
	683	#print "z",z
	684	#print "answer",answer
	685	assert num.allclose(z, answer)
	686
	687	def test_interpolate_point_outside_of_mesh(self):
	688	#Test linear interpolation of known values at vertices to
	689	#new points inside a triangle
	690
	691	a = [0.0, 0.0]
	692	b = [0.0, 5.0]
	693	c = [5.0, 0.0]
	694	d = [5.0, 5.0]
	695
	696	vertices = [a, b, c, d]
	697	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	698
	699	#Far away point
	700	d7 = [-1., -1.]
	701
	702	point_coords = [ d7]
	703	interp = Interpolate(vertices, triangles)
	704
	705	#Known values at vertices
	706	#Functions are x+y, x+2y, 2x+y, x-y-5
	707	f = [ [0., 0., 0., -5.], # (0,0)
	708	[5., 10., 5., -10.], # (0,5)
	709	[5., 5., 10.0, 0.], # (5,0)
	710	[10., 15., 15., -5.]] # (5,5)
	711
	712	z = interp.interpolate(f, point_coords) #, verbose=True)
	713	answer = num.array([ [NAN, NAN, NAN, NAN]]) # (-1,-1)
	714
	715	#print "***********"
	716	#print "z",z
	717	#print "answer",answer
	718	#print "***********"
	719
	720	#Should an error message be returned if points are outside
	721	# of the mesh?
	722	# A warning message is printed, if verbose is on.
	723
	724	for i in range(4):
	725	self.failUnless( z[0,i] == answer[0,i], 'Fail!')
	726
	727	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	728
	729	#print "z",z
	730	#print "answer",answer
	731
	732	for i in range(4):
	733	self.failUnless( z[0,i] == answer[0,i], 'Fail!')
	734
	735
	736	def test_interpolate_attributes_to_pointsIV(self):
	737	a = [-1.0, 0.0]
	738	b = [3.0, 4.0]
	739	c = [4.0, 1.0]
	740	d = [-3.0, 2.0] #3
	741	e = [-1.0, -2.0]
	742	f = [1.0, -2.0] #5
	743
	744	vertices = [a, b, c, d,e,f]
	745	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	746
	747
	748	point_coords = [[-2.0, 2.0],
	749	[-1.0, 1.0],
	750	[0.0, 2.0],
	751	[1.0, 1.0],
	752	[2.0, 1.0],
	753	[0.0, 0.0],
	754	[1.0, 0.0],
	755	[0.0, -1.0],
	756	[-0.2, -0.5],
	757	[-0.9, -1.5],
	758	[0.5, -1.9],
	759	[3.0, 1.0]]
	760
	761	interp = Interpolate(vertices, triangles)
	762	f = num.array([linear_function(vertices),2*linear_function(vertices)])
	763	f = num.transpose(f)
	764	#print "f",f
	765	z = interp.interpolate(f, point_coords)
	766	answer = [linear_function(point_coords),
	767	2*linear_function(point_coords) ]
	768	answer = num.transpose(answer)
	769	#print "z",z
	770	#print "answer",answer
	771	assert num.allclose(z, answer)
	772
	773	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	774
	775	#print "z",z
	776	#print "answer",answer
	777	assert num.allclose(z, answer)
	778
	779	def test_interpolate_blocking(self):
	780	a = [-1.0, 0.0]
	781	b = [3.0, 4.0]
	782	c = [4.0, 1.0]
	783	d = [-3.0, 2.0] #3
	784	e = [-1.0, -2.0]
	785	f = [1.0, -2.0] #5
	786
	787	vertices = [a, b, c, d,e,f]
	788	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	789
	790
	791	point_coords = [[-2.0, 2.0],
	792	[-1.0, 1.0],
	793	[0.0, 2.0],
	794	[1.0, 1.0],
	795	[2.0, 1.0],
	796	[0.0, 0.0],
	797	[1.0, 0.0],
	798	[0.0, -1.0],
	799	[-0.2, -0.5],
	800	[-0.9, -1.5],
	801	[0.5, -1.9],
	802	[3.0, 1.0]]
	803
	804	interp = Interpolate(vertices, triangles)
	805	f = num.array([linear_function(vertices),2*linear_function(vertices)])
	806	f = num.transpose(f)
	807	#print "f",f
	808	for blocking_max in range(len(point_coords)+2):
	809	#if True:
	810	# blocking_max = 5
	811	z = interp.interpolate(f, point_coords,
	812	start_blocking_len=blocking_max)
	813	answer = [linear_function(point_coords),
	814	2*linear_function(point_coords) ]
	815	answer = num.transpose(answer)
	816	#print "z",z
	817	#print "answer",answer
	818	assert num.allclose(z, answer)
	819
	820	f = num.array([linear_function(vertices),2*linear_function(vertices),
	821	2*linear_function(vertices) - 100])
	822	f = num.transpose(f)
	823	#print "f",f
	824	for blocking_max in range(len(point_coords)+2):
	825	#if True:
	826	# blocking_max = 5
	827	z = interp.interpolate(f, point_coords,
	828	start_blocking_len=blocking_max)
	829	answer = num.array([linear_function(point_coords),
	830	2*linear_function(point_coords) ,
	831	2*linear_function(point_coords)-100])
	832	z = num.transpose(z)
	833	#print "z",z
	834	#print "answer",answer
	835	assert num.allclose(z, answer)
	836
	837	def test_interpolate_geo_spatial(self):
	838	a = [-1.0, 0.0]
	839	b = [3.0, 4.0]
	840	c = [4.0, 1.0]
	841	d = [-3.0, 2.0] #3
	842	e = [-1.0, -2.0]
	843	f = [1.0, -2.0] #5
	844
	845	vertices = [a, b, c, d,e,f]
	846	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	847
	848
	849	point_coords_absolute = [[-2.0, 2.0],
	850	[-1.0, 1.0],
	851	[0.0, 2.0],
	852	[1.0, 1.0],
	853	[2.0, 1.0],
	854	[0.0, 0.0],
	855	[1.0, 0.0],
	856	[0.0, -1.0],
	857	[-0.2, -0.5],
	858	[-0.9, -1.5],
	859	[0.5, -1.9],
	860	[3.0, 1.0]]
	861
	862	geo = Geo_reference(57,100, 500)
	863	point_coords = geo.change_points_geo_ref(point_coords_absolute)
	864	point_coords = Geospatial_data(point_coords,geo_reference = geo)
	865
	866	interp = Interpolate(vertices, triangles)
	867	f = num.array([linear_function(vertices),2*linear_function(vertices)])
	868	f = num.transpose(f)
	869	#print "f",f
	870	for blocking_max in range(14):
	871	#if True:
	872	# blocking_max = 5
	873	z = interp.interpolate(f, point_coords,
	874	start_blocking_len=blocking_max)
	875	answer = [linear_function(point_coords.get_data_points( \
	876	absolute = True)),
	877	2*linear_function(point_coords.get_data_points( \
	878	absolute = True)) ]
	879	answer = num.transpose(answer)
	880	#print "z",z
	881	#print "answer",answer
	882	assert num.allclose(z, answer)
	883
	884	f = num.array([linear_function(vertices),2*linear_function(vertices),
	885	2*linear_function(vertices) - 100])
	886	f = num.transpose(f)
	887	#print "f",f
	888	for blocking_max in range(14):
	889	#if True:
	890	# blocking_max = 5
	891	z = interp.interpolate(f, point_coords,
	892	start_blocking_len=blocking_max)
	893	answer = num.array([linear_function(point_coords.get_data_points( \
	894	absolute = True)),
	895	2*linear_function(point_coords.get_data_points( \
	896	absolute = True)) ,
	897	2*linear_function(point_coords.get_data_points( \
	898	absolute = True))-100])
	899	z = num.transpose(z)
	900	#print "z",z
	901	#print "answer",answer
	902	assert num.allclose(z, answer)
	903
	904	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	905
	906	#print "z",z
	907	#print "answer",answer
	908	assert num.allclose(z, answer)
	909
	910	def test_interpolate_geo_spatial(self):
	911	a = [-1.0, 0.0]
	912	b = [3.0, 4.0]
	913	c = [4.0, 1.0]
	914	d = [-3.0, 2.0] #3
	915	e = [-1.0, -2.0]
	916	f = [1.0, -2.0] #5
	917
	918	vertices = [a, b, c, d,e,f]
	919	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	920
[6689]	921	point_coords_absolute = [[-2.0, 2.0],
	922	[-1.0, 1.0],
	923	[ 0.0, 2.0],
	924	[ 1.0, 1.0],
	925	[ 2.0, 1.0],
	926	[ 0.0, 0.0],
	927	[ 1.0, 0.0],
	928	[ 0.0, -1.0],
[6360]	929	[-0.2, -0.5],
	930	[-0.9, -1.5],
[6689]	931	[ 0.5, -1.9],
	932	[ 3.0, 1.0]]
[6174]	933
[6360]	934	geo = Geo_reference(57, 100, 500)
[6174]	935	point_coords = geo.change_points_geo_ref(point_coords_absolute)
[6360]	936	point_coords = Geospatial_data(point_coords, geo_reference=geo)
[6174]	937
	938	interp = Interpolate(vertices, triangles)
[6360]	939	f = num.array([linear_function(vertices), 2*linear_function(vertices)])
[6174]	940	f = num.transpose(f)
	941	z = interp.interpolate_block(f, point_coords)
[6360]	942	answer = [linear_function(point_coords.get_data_points(absolute=True)),
	943	2*linear_function(point_coords.get_data_points(absolute=True))
	944	]
[6174]	945	answer = num.transpose(answer)
[6360]	946	msg = ('Expected z\n%s\nto be close to answer\n%s'
	947	% (str(z), str(answer)))
	948	assert num.allclose(z, answer), msg
[6174]	949
	950	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	951
[6360]	952	msg = ('Expected z\n%s\nto be close to answer\n%s'
	953	% (str(z), str(answer)))
[6174]	954	assert num.allclose(z, answer)
	955
	956
	957	def test_interpolate_reuse_if_None(self):
	958	a = [-1.0, 0.0]
	959	b = [3.0, 4.0]
	960	c = [4.0, 1.0]
	961	d = [-3.0, 2.0] #3
	962	e = [-1.0, -2.0]
	963	f = [1.0, -2.0] #5
	964
	965	vertices = [a, b, c, d,e,f]
	966	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	967
	968
[6689]	969	point_coords = [[-2.0, 2.0],
	970	[-1.0, 1.0],
	971	[ 0.0, 2.0],
	972	[ 1.0, 1.0],
	973	[ 2.0, 1.0],
	974	[ 0.0, 0.0],
	975	[ 1.0, 0.0],
	976	[ 0.0, -1.0],
[6174]	977	[-0.2, -0.5],
	978	[-0.9, -1.5],
[6689]	979	[ 0.5, -1.9],
	980	[ 3.0, 1.0]]
[6174]	981
	982	interp = Interpolate(vertices, triangles)
	983	f = num.array([linear_function(vertices),2*linear_function(vertices)])
	984	f = num.transpose(f)
	985	z = interp.interpolate(f, point_coords,
	986	start_blocking_len=20)
	987	answer = [linear_function(point_coords),
	988	2*linear_function(point_coords) ]
	989	answer = num.transpose(answer)
	990	#print "z",z
	991	#print "answer",answer
	992	assert num.allclose(z, answer)
	993	assert num.allclose(interp._A_can_be_reused, True)
	994
	995	z = interp.interpolate(f)
	996	assert num.allclose(z, answer)
	997
	998	# This causes blocking to occur.
	999	z = interp.interpolate(f, start_blocking_len=10)
	1000	assert num.allclose(z, answer)
	1001	assert num.allclose(interp._A_can_be_reused, False)
	1002
	1003	#A is recalculated
	1004	z = interp.interpolate(f)
	1005	assert num.allclose(z, answer)
	1006	assert num.allclose(interp._A_can_be_reused, True)
	1007
	1008	interp = Interpolate(vertices, triangles)
	1009	#Must raise an exception, no points specified
	1010	try:
	1011	z = interp.interpolate(f)
	1012	except:
	1013	pass
	1014
	1015	def xxtest_interpolate_reuse_if_same(self):
	1016
	1017	# This on tests that repeated identical interpolation
	1018	# points makes use of precomputed matrix (Ole)
	1019	# This is not really a test and is disabled for now
	1020
	1021	a = [-1.0, 0.0]
	1022	b = [3.0, 4.0]
	1023	c = [4.0, 1.0]
	1024	d = [-3.0, 2.0] #3
	1025	e = [-1.0, -2.0]
	1026	f = [1.0, -2.0] #5
	1027
	1028	vertices = [a, b, c, d,e,f]
	1029	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	1030
	1031
[6689]	1032	point_coords = [[-2.0, 2.0],
	1033	[-1.0, 1.0],
	1034	[ 0.0, 2.0],
	1035	[ 1.0, 1.0],
	1036	[ 2.0, 1.0],
	1037	[ 0.0, 0.0],
	1038	[ 1.0, 0.0],
	1039	[ 0.0, -1.0],
[6174]	1040	[-0.2, -0.5],
	1041	[-0.9, -1.5],
[6689]	1042	[ 0.5, -1.9],
	1043	[ 3.0, 1.0]]
[6174]	1044
	1045	interp = Interpolate(vertices, triangles)
	1046	f = num.array([linear_function(vertices), 2*linear_function(vertices)])
	1047	f = num.transpose(f)
	1048	z = interp.interpolate(f, point_coords)
	1049	answer = [linear_function(point_coords),
	1050	2*linear_function(point_coords) ]
	1051	answer = num.transpose(answer)
	1052
	1053	assert num.allclose(z, answer)
	1054	assert num.allclose(interp._A_can_be_reused, True)
	1055
	1056
	1057	z = interp.interpolate(f) # None
	1058	assert num.allclose(z, answer)
	1059	z = interp.interpolate(f, point_coords) # Repeated (not really a test)
	1060	assert num.allclose(z, answer)
	1061
	1062
	1063
	1064	def test_interpolation_interface_time_only(self):
	1065
	1066	# Test spatio-temporal interpolation
	1067	# Test that spatio temporal function performs the correct
	1068	# interpolations in both time and space
	1069
	1070
	1071
	1072	#Three timesteps
	1073	time = [1.0, 5.0, 6.0]
	1074
	1075
	1076	#One quantity
[6304]	1077	Q = num.zeros( (3,6), num.float )
[6174]	1078
	1079	#Linear in time and space
	1080	a = [0.0, 0.0]
	1081	b = [0.0, 2.0]
	1082	c = [2.0, 0.0]
	1083	d = [0.0, 4.0]
	1084	e = [2.0, 2.0]
	1085	f = [4.0, 0.0]
	1086
	1087	points = [a, b, c, d, e, f]
	1088
	1089	for i, t in enumerate(time):
	1090	Q[i, :] = t*linear_function(points)
	1091
	1092
	1093	#Check basic interpolation of one quantity using averaging
	1094	#(no interpolation points or spatial info)
	1095	I = Interpolation_function(time, [mean(Q[0,:]),
	1096	mean(Q[1,:]),
	1097	mean(Q[2,:])])
	1098
	1099
	1100
	1101	#Check temporal interpolation
	1102	for i in [0,1,2]:
	1103	assert num.allclose(I(time[i]), mean(Q[i,:]))
	1104
	1105	#Midway
	1106	assert num.allclose(I( (time[0] + time[1])/2 ),
	1107	(I(time[0]) + I(time[1]))/2 )
	1108
	1109	assert num.allclose(I( (time[1] + time[2])/2 ),
	1110	(I(time[1]) + I(time[2]))/2 )
	1111
	1112	assert num.allclose(I( (time[0] + time[2])/2 ),
	1113	(I(time[0]) + I(time[2]))/2 )
	1114
	1115	#1/3
	1116	assert num.allclose(I( (time[0] + time[2])/3 ),
	1117	(I(time[0]) + I(time[2]))/3 )
	1118
	1119
	1120	#Out of bounds checks
	1121	try:
	1122	I(time[0]-1)
	1123	except:
	1124	pass
	1125	else:
	1126	raise 'Should raise exception'
	1127
	1128	try:
	1129	I(time[-1]+1)
	1130	except:
	1131	pass
	1132	else:
	1133	raise 'Should raise exception'
	1134
	1135
	1136
	1137
	1138	def test_interpolation_interface_spatial_only(self):
	1139	# Test spatio-temporal interpolation with constant time
	1140
	1141	#Three timesteps
	1142	time = [1.0, 5.0, 6.0]
	1143
	1144	#Setup mesh used to represent fitted function
	1145	a = [0.0, 0.0]
	1146	b = [0.0, 2.0]
	1147	c = [2.0, 0.0]
	1148	d = [0.0, 4.0]
	1149	e = [2.0, 2.0]
	1150	f = [4.0, 0.0]
	1151
	1152	points = [a, b, c, d, e, f]
	1153	#bac, bce, ecf, dbe
	1154	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1155
	1156
	1157	#New datapoints where interpolated values are sought
	1158	interpolation_points = [[ 0.0, 0.0],
	1159	[ 0.5, 0.5],
	1160	[ 0.7, 0.7],
	1161	[ 1.0, 0.5],
	1162	[ 2.0, 0.4],
	1163	[ 2.8, 1.2]]
	1164
	1165
	1166	#One quantity linear in space
	1167	Q = linear_function(points)
	1168
	1169
	1170	#Check interpolation of one quantity using interpolaton points
	1171	I = Interpolation_function(time, Q,
	1172	vertex_coordinates = points,
	1173	triangles = triangles,
	1174	interpolation_points = interpolation_points,
	1175	verbose = False)
	1176
	1177
	1178	answer = linear_function(interpolation_points)
	1179
	1180	t = time[0]
	1181	for j in range(50): #t in [1, 6]
	1182	for id in range(len(interpolation_points)):
	1183	assert num.allclose(I(t, id), answer[id])
	1184	t += 0.1
	1185
	1186	try:
	1187	I(1)
	1188	except:
	1189	pass
	1190	else:
	1191	raise 'Should raise exception'
	1192
	1193
	1194	def test_interpolation_interface(self):
	1195	# Test spatio-temporal interpolation
	1196	# Test that spatio temporal function performs the correct
	1197	# interpolations in both time and space
	1198
	1199	#Three timesteps
	1200	time = [1.0, 5.0, 6.0]
	1201
	1202	#Setup mesh used to represent fitted function
	1203	a = [0.0, 0.0]
	1204	b = [0.0, 2.0]
	1205	c = [2.0, 0.0]
	1206	d = [0.0, 4.0]
	1207	e = [2.0, 2.0]
	1208	f = [4.0, 0.0]
	1209
	1210	points = [a, b, c, d, e, f]
	1211	#bac, bce, ecf, dbe
	1212	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1213
	1214
	1215	#New datapoints where interpolated values are sought
	1216	interpolation_points = [[ 0.0, 0.0],
	1217	[ 0.5, 0.5],
	1218	[ 0.7, 0.7],
	1219	[ 1.0, 0.5],
	1220	[ 2.0, 0.4],
	1221	[ 2.8, 1.2]]
	1222
	1223	#One quantity
[6304]	1224	Q = num.zeros( (3,6), num.float )
[6174]	1225
	1226	#Linear in time and space
	1227	for i, t in enumerate(time):
	1228	Q[i, :] = t*linear_function(points)
	1229
	1230	#Check interpolation of one quantity using interpolaton points)
	1231	I = Interpolation_function(time, Q,
	1232	vertex_coordinates = points,
	1233	triangles = triangles,
	1234	interpolation_points = interpolation_points,
	1235	verbose = False)
	1236
	1237	answer = linear_function(interpolation_points)
	1238
	1239	t = time[0]
	1240	for j in range(50): #t in [1, 6]
	1241	for id in range(len(interpolation_points)):
	1242	assert num.allclose(I(t, id), t*answer[id])
	1243	t += 0.1
	1244
	1245	try:
	1246	I(1)
	1247	except:
	1248	pass
	1249	else:
	1250	raise 'Should raise exception'
	1251
	1252
	1253
	1254	def test_interpolation_interface_with_time_thinning(self):
	1255	# Test spatio-temporal interpolation
	1256	# Test that spatio temporal function performs the correct
	1257	# interpolations in both time and space
	1258
	1259	# Eight timesteps
	1260	time = [1.0, 2.0, 4.0, 5.0, 7.0, 8.0, 9.0, 10.0]
	1261
	1262	# Setup mesh used to represent fitted function
	1263	a = [0.0, 0.0]
	1264	b = [0.0, 2.0]
	1265	c = [2.0, 0.0]
	1266	d = [0.0, 4.0]
	1267	e = [2.0, 2.0]
	1268	f = [4.0, 0.0]
	1269
	1270	points = [a, b, c, d, e, f]
	1271	# bac, bce, ecf, dbe
	1272	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1273
	1274
	1275	# New datapoints where interpolated values are sought
	1276	interpolation_points = [[ 0.0, 0.0],
	1277	[ 0.5, 0.5],
	1278	[ 0.7, 0.7],
	1279	[ 1.0, 0.5],
	1280	[ 2.0, 0.4],
	1281	[ 2.8, 1.2]]
	1282
	1283	# One quantity
[6304]	1284	Q = num.zeros((8,6), num.float)
[6174]	1285
	1286	# Linear in time and space
	1287	for i, t in enumerate(time):
	1288	Q[i, :] = t*linear_function(points)
	1289
	1290	# Check interpolation of one quantity using interpolaton points) using default
	1291	# time_thinning of 1
	1292	I = Interpolation_function(time, Q,
	1293	vertex_coordinates=points,
	1294	triangles=triangles,
	1295	interpolation_points=interpolation_points,
	1296	verbose=False)
	1297
	1298	answer = linear_function(interpolation_points)
	1299
	1300
	1301	t = time[0]
	1302	for j in range(90): #t in [1, 10]
	1303	for id in range(len(interpolation_points)):
	1304	assert num.allclose(I(t, id), t*answer[id])
	1305	t += 0.1
	1306
	1307
	1308	# Now check time_thinning
	1309	I = Interpolation_function(time, Q,
	1310	vertex_coordinates=points,
	1311	triangles=triangles,
	1312	interpolation_points=interpolation_points,
	1313	time_thinning=2,
	1314	verbose=False)
	1315
	1316
	1317	assert len(I.time) == 4
	1318	assert( num.allclose(I.time, [1.0, 4.0, 7.0, 9.0] ))
	1319
	1320	answer = linear_function(interpolation_points)
	1321
	1322	t = time[0]
	1323	for j in range(80): #t in [1, 9]
	1324	for id in range(len(interpolation_points)):
	1325	assert num.allclose(I(t, id), t*answer[id])
	1326	t += 0.1
	1327
	1328
	1329
	1330	def test_interpolation_precompute_points(self):
	1331	# looking at a discrete mesh
	1332	#
	1333
	1334	#Three timesteps
	1335	time = [0.0, 60.0]
	1336
	1337	#Setup mesh used to represent fitted function
	1338	points = [[ 15., -20.],
	1339	[ 15., 10.],
	1340	[ 0., -20.],
	1341	[ 0., 10.],
	1342	[ 0., -20.],
	1343	[ 15., 10.]]
	1344
	1345	triangles = [[0, 1, 2],
	1346	[3, 4, 5]]
	1347
	1348	#New datapoints where interpolated values are sought
	1349	interpolation_points = [[ 1., 0.], [0.,1.]]
	1350
	1351	#One quantity
[6304]	1352	Q = num.zeros( (2,6), num.float )
[6174]	1353
	1354	#Linear in time and space
	1355	for i, t in enumerate(time):
	1356	Q[i, :] = t*linear_function(points)
	1357	#print "Q", Q
	1358
	1359
	1360
	1361	interp = Interpolate(points, triangles)
	1362	f = num.array([linear_function(points),2*linear_function(points)])
	1363	f = num.transpose(f)
	1364	#print "f",f
	1365	z = interp.interpolate(f, interpolation_points)
	1366	answer = [linear_function(interpolation_points),
	1367	2*linear_function(interpolation_points) ]
	1368	answer = num.transpose(answer)
	1369	#print "z",z
	1370	#print "answer",answer
	1371	assert num.allclose(z, answer)
	1372
	1373
	1374	#Check interpolation of one quantity using interpolaton points)
	1375	I = Interpolation_function(time, Q,
	1376	vertex_coordinates = points,
	1377	triangles = triangles,
	1378	interpolation_points = interpolation_points,
	1379	verbose = False)
	1380
	1381	#print "I.precomputed_values", I.precomputed_values
	1382
	1383	msg = 'Interpolation failed'
	1384	assert num.allclose(I.precomputed_values['Attribute'][1], [60, 60]), msg
	1385	#self.failUnless( I.precomputed_values['Attribute'][1] == 60.0,
	1386	# ' failed')
	1387
	1388	def test_interpolation_function_outside_point(self):
	1389	# Test spatio-temporal interpolation
	1390	# Test that spatio temporal function performs the correct
	1391	# interpolations in both time and space
	1392
	1393	# Three timesteps
	1394	time = [1.0, 5.0, 6.0]
	1395
	1396	# Setup mesh used to represent fitted function
	1397	a = [0.0, 0.0]
	1398	b = [0.0, 2.0]
	1399	c = [2.0, 0.0]
	1400	d = [0.0, 4.0]
	1401	e = [2.0, 2.0]
	1402	f = [4.0, 0.0]
	1403
	1404	points = [a, b, c, d, e, f]
	1405	#bac, bce, ecf, dbe
	1406	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1407
	1408
	1409	# New datapoints where interpolated values are sought
	1410	interpolation_points = [[ 0.0, 0.0],
	1411	[ 0.5, 0.5],
	1412	[ 0.7, 0.7],
	1413	[ 1.0, 0.5],
	1414	[ 2.0, 0.4],
	1415	[ 545354534, 4354354353]] # outside the mesh
	1416
	1417	# One quantity
[6304]	1418	Q = num.zeros( (3,6), num.float )
[6174]	1419
	1420	# Linear in time and space
	1421	for i, t in enumerate(time):
	1422	Q[i, :] = t*linear_function(points)
	1423
	1424	# Check interpolation of one quantity using interpolaton points)
	1425
	1426	I = Interpolation_function(time, Q,
	1427	vertex_coordinates = points,
	1428	triangles = triangles,
	1429	interpolation_points = interpolation_points,
	1430	verbose = False)
	1431
	1432
	1433	assert num.alltrue(I.precomputed_values['Attribute'][:,4] != NAN)
	1434	assert num.sometrue(I.precomputed_values['Attribute'][:,5] == NAN)
	1435
	1436	#X = I.precomputed_values['Attribute'][1,:]
	1437	#print X
	1438	#print take(X, X == NAN)
	1439	#print where(X == NAN, range(len(X)), 0)
	1440
	1441	answer = linear_function(interpolation_points)
	1442
	1443	t = time[0]
	1444	for j in range(50): #t in [1, 6]
	1445	for id in range(len(interpolation_points)-1):
	1446	assert num.allclose(I(t, id), t*answer[id])
	1447	t += 0.1
	1448
	1449	# Now test the point outside the mesh
	1450	t = time[0]
	1451	for j in range(50): #t in [1, 6]
	1452	self.failUnless(I(t, 5) == NAN, 'Fail!')
	1453	t += 0.1
	1454
	1455	try:
	1456	I(1)
	1457	except:
	1458	pass
	1459	else:
	1460	raise 'Should raise exception'
	1461
	1462
	1463	def test_interpolation_function_time(self):
	1464	#Test a long time series with an error in it (this did cause an
	1465	#error once)
	1466
	1467
	1468	time = num.array(\
	1469	[0.00000000e+00, 5.00000000e-02, 1.00000000e-01, 1.50000000e-01,
	1470	2.00000000e-01, 2.50000000e-01, 3.00000000e-01, 3.50000000e-01,
	1471	4.00000000e-01, 4.50000000e-01, 5.00000000e-01, 5.50000000e-01,
	1472	6.00000000e-01, 6.50000000e-01, 7.00000000e-01, 7.50000000e-01,
	1473	8.00000000e-01, 8.50000000e-01, 9.00000000e-01, 9.50000000e-01,
	1474	1.00000000e-00, 1.05000000e+00, 1.10000000e+00, 1.15000000e+00,
	1475	1.20000000e+00, 1.25000000e+00, 1.30000000e+00, 1.35000000e+00,
	1476	1.40000000e+00, 1.45000000e+00, 1.50000000e+00, 1.55000000e+00,
	1477	1.60000000e+00, 1.65000000e+00, 1.70000000e+00, 1.75000000e+00,
	1478	1.80000000e+00, 1.85000000e+00, 1.90000000e+00, 1.95000000e+00,
	1479	2.00000000e+00, 2.05000000e+00, 2.10000000e+00, 2.15000000e+00,
	1480	2.20000000e+00, 2.25000000e+00, 2.30000000e+00, 2.35000000e+00,
	1481	2.40000000e+00, 2.45000000e+00, 2.50000000e+00, 2.55000000e+00,
	1482	2.60000000e+00, 2.65000000e+00, 2.70000000e+00, 2.75000000e+00,
	1483	2.80000000e+00, 2.85000000e+00, 2.90000000e+00, 2.95000000e+00,
	1484	3.00000000e+00, 3.05000000e+00, 9.96920997e+36, 3.15000000e+00,
	1485	3.20000000e+00, 3.25000000e+00, 3.30000000e+00, 3.35000000e+00,
	1486	3.40000000e+00, 3.45000000e+00, 3.50000000e+00, 3.55000000e+00,
	1487	3.60000000e+00, 3.65000000e+00, 3.70000000e+00, 3.75000000e+00,
	1488	3.80000000e+00, 3.85000000e+00, 3.90000000e+00, 3.95000000e+00,
	1489	4.00000000e+00, 4.05000000e+00, 4.10000000e+00, 4.15000000e+00,
	1490	4.20000000e+00, 4.25000000e+00, 4.30000000e+00, 4.35000000e+00,
	1491	4.40000000e+00, 4.45000000e+00, 4.50000000e+00, 4.55000000e+00,
	1492	4.60000000e+00, 4.65000000e+00, 4.70000000e+00, 4.75000000e+00,
	1493	4.80000000e+00, 4.85000000e+00, 4.90000000e+00, 4.95000000e+00,
	1494	5.00000000e+00, 5.05000000e+00, 5.10000000e+00, 5.15000000e+00,
	1495	5.20000000e+00, 5.25000000e+00, 5.30000000e+00, 5.35000000e+00,
	1496	5.40000000e+00, 5.45000000e+00, 5.50000000e+00, 5.55000000e+00,
	1497	5.60000000e+00, 5.65000000e+00, 5.70000000e+00, 5.75000000e+00,
	1498	5.80000000e+00, 5.85000000e+00, 5.90000000e+00, 5.95000000e+00,
	1499	6.00000000e+00, 6.05000000e+00, 6.10000000e+00, 6.15000000e+00,
	1500	6.20000000e+00, 6.25000000e+00, 6.30000000e+00, 6.35000000e+00,
	1501	6.40000000e+00, 6.45000000e+00, 6.50000000e+00, 6.55000000e+00,
	1502	6.60000000e+00, 6.65000000e+00, 6.70000000e+00, 6.75000000e+00,
	1503	6.80000000e+00, 6.85000000e+00, 6.90000000e+00, 6.95000000e+00,
	1504	7.00000000e+00, 7.05000000e+00, 7.10000000e+00, 7.15000000e+00,
	1505	7.20000000e+00, 7.25000000e+00, 7.30000000e+00, 7.35000000e+00,
	1506	7.40000000e+00, 7.45000000e+00, 7.50000000e+00, 7.55000000e+00,
	1507	7.60000000e+00, 7.65000000e+00, 7.70000000e+00, 7.75000000e+00,
	1508	7.80000000e+00, 7.85000000e+00, 7.90000000e+00, 7.95000000e+00,
	1509	8.00000000e+00, 8.05000000e+00, 8.10000000e+00, 8.15000000e+00,
	1510	8.20000000e+00, 8.25000000e+00, 8.30000000e+00, 8.35000000e+00,
	1511	8.40000000e+00, 8.45000000e+00, 8.50000000e+00, 8.55000000e+00,
	1512	8.60000000e+00, 8.65000000e+00, 8.70000000e+00, 8.75000000e+00,
	1513	8.80000000e+00, 8.85000000e+00, 8.90000000e+00, 8.95000000e+00,
	1514	9.00000000e+00, 9.05000000e+00, 9.10000000e+00, 9.15000000e+00,
	1515	9.20000000e+00, 9.25000000e+00, 9.30000000e+00, 9.35000000e+00,
	1516	9.40000000e+00, 9.45000000e+00, 9.50000000e+00, 9.55000000e+00,
	1517	9.60000000e+00, 9.65000000e+00, 9.70000000e+00, 9.75000000e+00,
	1518	9.80000000e+00, 9.85000000e+00, 9.90000000e+00, 9.95000000e+00,
	1519	1.00000000e+01, 1.00500000e+01, 1.01000000e+01, 1.01500000e+01,
	1520	1.02000000e+01, 1.02500000e+01, 1.03000000e+01, 1.03500000e+01,
	1521	1.04000000e+01, 1.04500000e+01, 1.05000000e+01, 1.05500000e+01,
	1522	1.06000000e+01, 1.06500000e+01, 1.07000000e+01, 1.07500000e+01,
	1523	1.08000000e+01, 1.08500000e+01, 1.09000000e+01, 1.09500000e+01,
	1524	1.10000000e+01, 1.10500000e+01, 1.11000000e+01, 1.11500000e+01,
	1525	1.12000000e+01, 1.12500000e+01, 1.13000000e+01, 1.13500000e+01,
	1526	1.14000000e+01, 1.14500000e+01, 1.15000000e+01, 1.15500000e+01,
	1527	1.16000000e+01, 1.16500000e+01, 1.17000000e+01, 1.17500000e+01,
	1528	1.18000000e+01, 1.18500000e+01, 1.19000000e+01, 1.19500000e+01,
	1529	1.20000000e+01, 1.20500000e+01, 1.21000000e+01, 1.21500000e+01,
	1530	1.22000000e+01, 1.22500000e+01, 1.23000000e+01, 1.23500000e+01,
	1531	1.24000000e+01, 1.24500000e+01, 1.25000000e+01, 1.25500000e+01,
	1532	1.26000000e+01, 1.26500000e+01, 1.27000000e+01, 1.27500000e+01,
	1533	1.28000000e+01, 1.28500000e+01, 1.29000000e+01, 1.29500000e+01,
	1534	1.30000000e+01, 1.30500000e+01, 1.31000000e+01, 1.31500000e+01,
	1535	1.32000000e+01, 1.32500000e+01, 1.33000000e+01, 1.33500000e+01,
	1536	1.34000000e+01, 1.34500000e+01, 1.35000000e+01, 1.35500000e+01,
	1537	1.36000000e+01, 1.36500000e+01, 1.37000000e+01, 1.37500000e+01,
	1538	1.38000000e+01, 1.38500000e+01, 1.39000000e+01, 1.39500000e+01,
	1539	1.40000000e+01, 1.40500000e+01, 1.41000000e+01, 1.41500000e+01,
	1540	1.42000000e+01, 1.42500000e+01, 1.43000000e+01, 1.43500000e+01,
	1541	1.44000000e+01, 1.44500000e+01, 1.45000000e+01, 1.45500000e+01,
	1542	1.46000000e+01, 1.46500000e+01, 1.47000000e+01, 1.47500000e+01,
	1543	1.48000000e+01, 1.48500000e+01, 1.49000000e+01, 1.49500000e+01,
	1544	1.50000000e+01, 1.50500000e+01, 1.51000000e+01, 1.51500000e+01,
	1545	1.52000000e+01, 1.52500000e+01, 1.53000000e+01, 1.53500000e+01,
	1546	1.54000000e+01, 1.54500000e+01, 1.55000000e+01, 1.55500000e+01,
	1547	1.56000000e+01, 1.56500000e+01, 1.57000000e+01, 1.57500000e+01,
	1548	1.58000000e+01, 1.58500000e+01, 1.59000000e+01, 1.59500000e+01,
	1549	1.60000000e+01, 1.60500000e+01, 1.61000000e+01, 1.61500000e+01,
	1550	1.62000000e+01, 1.62500000e+01, 1.63000000e+01, 1.63500000e+01,
	1551	1.64000000e+01, 1.64500000e+01, 1.65000000e+01, 1.65500000e+01,
	1552	1.66000000e+01, 1.66500000e+01, 1.67000000e+01, 1.67500000e+01,
	1553	1.68000000e+01, 1.68500000e+01, 1.69000000e+01, 1.69500000e+01,
	1554	1.70000000e+01, 1.70500000e+01, 1.71000000e+01, 1.71500000e+01,
	1555	1.72000000e+01, 1.72500000e+01, 1.73000000e+01, 1.73500000e+01,
	1556	1.74000000e+01, 1.74500000e+01, 1.75000000e+01, 1.75500000e+01,
	1557	1.76000000e+01, 1.76500000e+01, 1.77000000e+01, 1.77500000e+01,
	1558	1.78000000e+01, 1.78500000e+01, 1.79000000e+01, 1.79500000e+01,
	1559	1.80000000e+01, 1.80500000e+01, 1.81000000e+01, 1.81500000e+01,
	1560	1.82000000e+01, 1.82500000e+01, 1.83000000e+01, 1.83500000e+01,
	1561	1.84000000e+01, 1.84500000e+01, 1.85000000e+01, 1.85500000e+01,
	1562	1.86000000e+01, 1.86500000e+01, 1.87000000e+01, 1.87500000e+01,
	1563	1.88000000e+01, 1.88500000e+01, 1.89000000e+01, 1.89500000e+01,
	1564	1.90000000e+01, 1.90500000e+01, 1.91000000e+01, 1.91500000e+01,
	1565	1.92000000e+01, 1.92500000e+01, 1.93000000e+01, 1.93500000e+01,
	1566	1.94000000e+01, 1.94500000e+01, 1.95000000e+01, 1.95500000e+01,
	1567	1.96000000e+01, 1.96500000e+01, 1.97000000e+01, 1.97500000e+01,
	1568	1.98000000e+01, 1.98500000e+01, 1.99000000e+01, 1.99500000e+01,
	1569	2.00000000e+01, 2.00500000e+01, 2.01000000e+01, 2.01500000e+01,
	1570	2.02000000e+01, 2.02500000e+01, 2.03000000e+01, 2.03500000e+01,
	1571	2.04000000e+01, 2.04500000e+01, 2.05000000e+01, 2.05500000e+01,
	1572	2.06000000e+01, 2.06500000e+01, 2.07000000e+01, 2.07500000e+01,
	1573	2.08000000e+01, 2.08500000e+01, 2.09000000e+01, 2.09500000e+01,
	1574	2.10000000e+01, 2.10500000e+01, 2.11000000e+01, 2.11500000e+01,
	1575	2.12000000e+01, 2.12500000e+01, 2.13000000e+01, 2.13500000e+01,
	1576	2.14000000e+01, 2.14500000e+01, 2.15000000e+01, 2.15500000e+01,
	1577	2.16000000e+01, 2.16500000e+01, 2.17000000e+01, 2.17500000e+01,
	1578	2.18000000e+01, 2.18500000e+01, 2.19000000e+01, 2.19500000e+01,
	1579	2.20000000e+01, 2.20500000e+01, 2.21000000e+01, 2.21500000e+01,
	1580	2.22000000e+01, 2.22500000e+01, 2.23000000e+01, 2.23500000e+01,
	1581	2.24000000e+01, 2.24500000e+01, 2.25000000e+01])
	1582
	1583	#print 'Diff', time[1:] - time[:-1]
	1584
	1585	#Setup mesh used to represent fitted function
	1586	a = [0.0, 0.0]
	1587	b = [0.0, 2.0]
	1588	c = [2.0, 0.0]
	1589	d = [0.0, 4.0]
	1590	e = [2.0, 2.0]
	1591	f = [4.0, 0.0]
	1592
	1593	points = [a, b, c, d, e, f]
	1594	#bac, bce, ecf, dbe
	1595	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1596
	1597
	1598	#New datapoints where interpolated values are sought
	1599	interpolation_points = [[ 0.0, 0.0],
	1600	[ 0.5, 0.5],
	1601	[ 0.7, 0.7],
	1602	[ 1.0, 0.5],
	1603	[ 2.0, 0.4],
	1604	[ 545354534, 4354354353]] # outside the mesh
	1605
	1606	#One quantity
[6304]	1607	Q = num.zeros( (len(time),6), num.float )
[6174]	1608
	1609	#Linear in time and space
	1610	for i, t in enumerate(time):
	1611	Q[i, :] = t*linear_function(points)
	1612
	1613	#Check interpolation of one quantity using interpolaton points)
	1614	try:
	1615	I = Interpolation_function(time, Q,
	1616	vertex_coordinates = points,
	1617	triangles = triangles,
	1618	interpolation_points = interpolation_points,
	1619	verbose = False)
	1620	except:
	1621	pass
	1622	else:
	1623	raise 'Should raise exception due to time being non-monotoneous'
	1624
	1625
	1626	def test_points_outside_the_polygon(self):
	1627	a = [-1.0, 0.0]
	1628	b = [3.0, 4.0]
	1629	c = [4.0, 1.0]
	1630	d = [-3.0, 2.0] #3
	1631	e = [-1.0, -2.0]
	1632	f = [1.0, -2.0] #5
	1633
	1634	vertices = [a, b, c, d,e,f]
	1635	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	1636
	1637	point_coords = [[-2.0, 2.0],
	1638	[-1.0, 1.0],
	1639	[9999.0, 9999.0], # point Outside poly
	1640	[-9999.0, 1.0], # point Outside poly
	1641	[2.0, 1.0],
	1642	[0.0, 0.0],
	1643	[1.0, 0.0],
	1644	[0.0, -1.0],
	1645	[-0.2, -0.5],
	1646	[-0.9, -1.5],
	1647	[0.5, -1.9],
	1648	[999999, 9999999]] # point Outside poly
	1649	geo_data = Geospatial_data(data_points = point_coords)
	1650
	1651	interp = Interpolate(vertices, triangles)
	1652	f = num.array([linear_function(vertices),2*linear_function(vertices)])
	1653	f = num.transpose(f)
	1654	#print "f",f
	1655	z = interp.interpolate(f, geo_data)
	1656	#z = interp.interpolate(f, point_coords)
	1657	answer = [linear_function(point_coords),
	1658	2*linear_function(point_coords) ]
	1659	answer = num.transpose(answer)
	1660	answer[2,:] = [NAN, NAN]
	1661	answer[3,:] = [NAN, NAN]
	1662	answer[11,:] = [NAN, NAN]
	1663	#print "z",z
	1664	#print "answer _ fixed",answer
	1665	assert num.allclose(z[0:1], answer[0:1])
	1666	assert num.allclose(z[4:10], answer[4:10])
	1667	for i in [2,3,11]:
	1668	self.failUnless( z[i,1] == answer[11,1], 'Fail!')
	1669	self.failUnless( z[i,0] == answer[11,0], 'Fail!')
	1670
	1671	def test_interpolate_sww2csv(self):
	1672
	1673	def elevation_function(x, y):
	1674	return -x
	1675
	1676	# Create mesh
	1677	mesh_file = tempfile.mktemp(".tsh")
	1678	points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
	1679	m = Mesh()
	1680	m.add_vertices(points)
	1681	m.auto_segment()
	1682	m.generate_mesh(verbose=False)
	1683	m.export_mesh_file(mesh_file)
	1684
	1685	# Create shallow water domain
	1686	domain = Domain(mesh_file)
	1687	os.remove(mesh_file)
	1688
	1689	domain.default_order = 2
	1690
	1691	# This test was made before tight_slope_limiters were introduced
	1692	# Since were are testing interpolation values this is OK
	1693	domain.tight_slope_limiters = 0
	1694
	1695	# Set some field values
	1696	domain.set_quantity('elevation', elevation_function)
	1697	domain.set_quantity('friction', 0.03)
	1698	domain.set_quantity('xmomentum', 3.0)
	1699	domain.set_quantity('ymomentum', 4.0)
	1700
	1701	######################
	1702	# Boundary conditions
	1703	B = Transmissive_boundary(domain)
	1704	domain.set_boundary( {'exterior': B})
	1705
	1706	# This call mangles the stage values.
	1707	domain.distribute_to_vertices_and_edges()
	1708	domain.set_quantity('stage', 1.0)
	1709
	1710
	1711	domain.set_name('datatest' + str(time.time()))
	1712	domain.format = 'sww'
	1713	domain.smooth = True
	1714	domain.reduction = mean
	1715
	1716	sww = get_dataobject(domain)
	1717	sww.store_connectivity()
	1718	sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
	1719	domain.set_quantity('stage', 10.0) # This is automatically limited
	1720	# So it will not be less than the elevation
	1721	domain.time = 2.
	1722	sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
	1723
	1724	# Test the function
	1725	points = [[5.0,1.],[0.5,2.]]
	1726	depth_file = tempfile.mktemp(".csv")
	1727	velocity_x_file = tempfile.mktemp(".csv")
	1728	velocity_y_file = tempfile.mktemp(".csv")
	1729	interpolate_sww2csv(sww.filename, points, depth_file,
	1730	velocity_x_file,
	1731	velocity_y_file,
	1732	verbose=False)
	1733
	1734	depth_answers_array = [[0.0, 6.0, 1.5], [2.0, 15., 10.5]]
	1735	velocity_x_answers_array = [[0.0, 3./6.0, 3./1.5],
	1736	[2.0, 3./15., 3/10.5]]
	1737	velocity_y_answers_array = [[0.0, 4./6.0, 4./1.5],
	1738	[2.0, 4./15., 4./10.5]]
	1739	depth_file_handle = file(depth_file)
	1740	depth_reader = csv.reader(depth_file_handle)
	1741	depth_reader.next()
	1742	velocity_x_file_handle = file(velocity_x_file)
	1743	velocity_x_reader = csv.reader(velocity_x_file_handle)
	1744	velocity_x_reader.next()
	1745	for depths, velocitys, depth_answers, velocity_answers in map(None,
	1746	depth_reader,
	1747	velocity_x_reader,
	1748	depth_answers_array,
	1749	velocity_x_answers_array):
	1750	for i in range(len(depths)):
	1751	#print "depths",depths[i]
	1752	#print "depth_answers",depth_answers[i]
	1753	#print "velocitys",velocitys[i]
	1754	#print "velocity_answers_array", velocity_answers[i]
	1755	msg = 'Interpolation failed'
	1756	assert num.allclose(float(depths[i]), depth_answers[i]), msg
	1757	assert num.allclose(float(velocitys[i]), velocity_answers[i]), msg
	1758
	1759	velocity_y_file_handle = file(velocity_y_file)
	1760	velocity_y_reader = csv.reader(velocity_y_file_handle)
	1761	velocity_y_reader.next()
	1762	for velocitys, velocity_answers in map(None,
	1763	velocity_y_reader,
	1764	velocity_y_answers_array):
	1765	for i in range(len(depths)):
	1766	#print "depths",depths[i]
	1767	#print "depth_answers",depth_answers[i]
	1768	#print "velocitys",velocitys[i]
	1769	#print "velocity_answers_array", velocity_answers[i]
	1770	msg = 'Interpolation failed'
	1771	assert num.allclose(float(depths[i]), depth_answers[i]), msg
	1772	assert num.allclose(float(velocitys[i]), velocity_answers[i]), msg
	1773
	1774	# clean up
	1775	depth_file_handle.close()
	1776	velocity_y_file_handle.close()
	1777	velocity_x_file_handle.close()
	1778	#print "sww.filename",sww.filename
	1779	os.remove(sww.filename)
	1780	os.remove(depth_file)
	1781	os.remove(velocity_x_file)
	1782	os.remove(velocity_y_file)
	1783
	1784
	1785	def test_interpolate_one_point_many_triangles(self):
	1786	# this test has 10 triangles that share the same vert.
	1787	# If the number of points per cell in a quad tree is less
	1788	# than 10 it will crash.
	1789	z0 = [2.0, 5.0]
	1790	z1 = [2.0, 5.0]
	1791	z2 = [2.0, 5.0]
	1792	z3 = [2.0, 5.0]
	1793	z4 = [2.0, 5.0]
	1794	z5 = [2.0, 5.0]
	1795	z6 = [2.0, 5.0]
	1796	z7 = [2.0, 5.0]
	1797	z8 = [2.0, 5.0]
	1798	z9 = [2.0, 5.0]
	1799	z10 = [2.0, 5.0]
	1800
	1801	v0 = [0.0, 0.0]
	1802	v1 = [1.0, 0.0]
	1803	v2 = [2.0, 0.0]
	1804	v3 = [3.0, 0.0]
	1805	v4 = [4.0, 0.0]
	1806	v5 = [0.0, 10.0]
	1807	v6 = [1.0, 10.0]
	1808	v7 = [2.0, 10.0]
	1809	v8 = [3.0, 10.0]
	1810	v9 = [4.0, 10.0]
	1811
	1812	vertices = [z0,v0, v1, v2, v3,v4 ,v5, v6, v7, v8, v9,
	1813	z1, z2, z3, z4, z5, z6, z7, z8, z9]
	1814	triangles = [
	1815	[11,1,2],
	1816	[12,2,3],
	1817	[13,3,4],
	1818	[14,4,5],
	1819	[7,6,15],
	1820	[8,7,16],
	1821	[9,8,17],
	1822	[10,9,18],
	1823	[6,1,19],
	1824	[5,10,0]
	1825	]
	1826
	1827	d0 = [1.0, 1.0]
	1828	d1 = [1.0, 2.0]
	1829	d2 = [3.0, 1.0]
	1830	point_coords = [ d0, d1, d2]
	1831	try:
	1832	interp = Interpolate(vertices, triangles)
	1833	except RuntimeError:
	1834	self.failUnless(0 ==1, 'quad fails with 10 verts at the same \
	1835	position. Real problems have had 9. \
	1836	Should be able to handle 13.')
	1837	f = linear_function(vertices)
	1838	z = interp.interpolate(f, point_coords)
	1839	answer = linear_function(point_coords)
	1840
	1841	#print "z",z
	1842	#print "answer",answer
	1843	assert num.allclose(z, answer)
[6360]	1844
	1845	################################################################################
	1846
[6174]	1847	if __name__ == "__main__":
	1848	suite = unittest.makeSuite(Test_Interpolate,'test')
[6304]	1849	runner = unittest.TextTestRunner() #verbosity=1)
[6174]	1850	runner.run(suite)
	1851

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

Context Navigation

source: branches/numpy/anuga/fit_interpolate/test_interpolate.py @ 6982

Download in other formats: