Context Navigation

source: anuga_work/development/anuga_1d/test_quantity.py @ 5858

Last change on this file since 5858 was 5858, checked in by steve, 16 years ago

File size: 82.5 KB

Rev	Line
[5536]	1	#!/usr/bin/env python
	2
	3	import unittest
	4	from math import sqrt, pi
	5	import tempfile
	6
	7	from quantity import *
[5858]	8	from config import epsilon
[5536]	9	from Numeric import allclose, array, ones, Float
[5858]	10
[5536]	11	from domain import Domain
	12
[5858]	13
[5536]	14	#Aux for fit_interpolate.fit example
	15	def linear_function(point):
	16	point = array(point)
	17	return point[:,0]+point[:,1]
	18
	19
	20	class Test_Quantity(unittest.TestCase):
	21	def setUp(self):
	22	from domain import Domain
	23
	24	a = [0.0, 0.0]
	25	b = [0.0, 2.0]
	26	c = [2.0, 0.0]
	27	d = [0.0, 4.0]
	28	e = [2.0, 2.0]
	29	f = [4.0, 0.0]
	30
	31	points = [a, b, c, d, e, f]
	32
	33	#bac, bce, ecf, dbe
	34	elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
	35
	36	self.mesh1 = Domain(points[:3], [elements[0]])
	37	self.mesh1.check_integrity()
	38
	39	self.mesh4 = Domain(points, elements)
	40	self.mesh4.check_integrity()
	41
	42	# UTM round Onslow
	43	a = [240000, 7620000]
	44	b = [240000, 7680000]
	45	c = [300000, 7620000]
	46
	47	points = [a, b, c]
	48	elements = [[0,2,1]]
	49
	50	self.mesh_onslow = Domain(points, elements)
	51	self.mesh_onslow.check_integrity()
	52
	53	def tearDown(self):
	54	pass
	55	#print " Tearing down"
	56
	57
	58	def test_creation(self):
	59
	60	quantity = Quantity(self.mesh1, [[1,2,3]])
	61	assert allclose(quantity.vertex_values, [[1.,2.,3.]])
	62
	63	try:
	64	quantity = Quantity()
	65	except:
	66	pass
	67	else:
	68	raise 'Should have raised empty quantity exception'
	69
	70
	71	try:
	72	quantity = Quantity([1,2,3])
	73	except AssertionError:
	74	pass
	75	except:
	76	raise 'Should have raised "mising mesh object" error'
	77
	78
	79	def test_creation_zeros(self):
	80
	81	quantity = Quantity(self.mesh1)
	82	assert allclose(quantity.vertex_values, [[0.,0.,0.]])
	83
	84
	85	quantity = Quantity(self.mesh4)
	86	assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
	87	[0.,0.,0.], [0.,0.,0.]])
	88
	89
	90	def test_interpolation(self):
	91	quantity = Quantity(self.mesh1, [[1,2,3]])
	92	assert allclose(quantity.centroid_values, [2.0]) #Centroid
	93
	94	assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
	95
	96
	97	def test_interpolation2(self):
	98	quantity = Quantity(self.mesh4,
	99	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	100	assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
	101
	102
	103	quantity.extrapolate_second_order()
	104
	105	#print quantity.vertex_values
	106	assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5],
	107	[3.+2./3, 6.+2./3, 4.+2./3],
	108	[4.6, 3.4, 1.],
	109	[-5.0, 1.0, 4.0]])
	110
	111	#print quantity.edge_values
	112	assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25],
	113	[5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6],
	114	[2.2, 2.8, 4.0],
	115	[2.5, -0.5, -2.0]])
	116
	117
	118	#assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
	119	# [5., 5., 5.],
	120	# [4.5, 4.5, 0.],
	121	# [3.0, -1.5, -1.5]])
	122
	123	def test_get_extrema_1(self):
	124	quantity = Quantity(self.mesh4,
	125	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	126	assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
	127
	128	v = quantity.get_maximum_value()
	129	assert v == 5
	130
	131	v = quantity.get_minimum_value()
	132	assert v == 0
	133
	134	i = quantity.get_maximum_index()
	135	assert i == 1
	136
	137	i = quantity.get_minimum_index()
	138	assert i == 3
	139
	140	x,y = quantity.get_maximum_location()
	141	xref, yref = 4.0/3, 4.0/3
	142	assert x == xref
	143	assert y == yref
	144
	145	v = quantity.get_values(interpolation_points = [[x,y]])
	146	assert allclose(v, 5)
	147
	148
	149	x,y = quantity.get_minimum_location()
	150	v = quantity.get_values(interpolation_points = [[x,y]])
	151	assert allclose(v, 0)
	152
	153
	154	def test_get_maximum_2(self):
	155
	156	a = [0.0, 0.0]
	157	b = [0.0, 2.0]
	158	c = [2.0,0.0]
	159	d = [0.0, 4.0]
	160	e = [2.0, 2.0]
	161	f = [4.0,0.0]
	162
	163	points = [a, b, c, d, e, f]
	164	#bac, bce, ecf, dbe
	165	vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	166
	167	domain = Domain(points, vertices)
	168
	169	quantity = Quantity(domain)
	170	quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
	171
	172	v = quantity.get_maximum_value()
	173	assert v == 6
	174
	175	v = quantity.get_minimum_value()
	176	assert v == 2
	177
	178	i = quantity.get_maximum_index()
	179	assert i == 3
	180
	181	i = quantity.get_minimum_index()
	182	assert i == 0
	183
	184	x,y = quantity.get_maximum_location()
	185	xref, yref = 2.0/3, 8.0/3
	186	assert x == xref
	187	assert y == yref
	188
	189	v = quantity.get_values(interpolation_points = [[x,y]])
	190	assert allclose(v, 6)
	191
	192	x,y = quantity.get_minimum_location()
	193	v = quantity.get_values(interpolation_points = [[x,y]])
	194	assert allclose(v, 2)
	195
	196	#Multiple locations for maximum -
	197	#Test that the algorithm picks the first occurrence
	198	v = quantity.get_maximum_value(indices=[0,1,2])
	199	assert allclose(v, 4)
	200
	201	i = quantity.get_maximum_index(indices=[0,1,2])
	202	assert i == 1
	203
	204	x,y = quantity.get_maximum_location(indices=[0,1,2])
	205	xref, yref = 4.0/3, 4.0/3
	206	assert x == xref
	207	assert y == yref
	208
	209	v = quantity.get_values(interpolation_points = [[x,y]])
	210	assert allclose(v, 4)
	211
	212	# More test of indices......
	213	v = quantity.get_maximum_value(indices=[2,3])
	214	assert allclose(v, 6)
	215
	216	i = quantity.get_maximum_index(indices=[2,3])
	217	assert i == 3
	218
	219	x,y = quantity.get_maximum_location(indices=[2,3])
	220	xref, yref = 2.0/3, 8.0/3
	221	assert x == xref
	222	assert y == yref
	223
	224	v = quantity.get_values(interpolation_points = [[x,y]])
	225	assert allclose(v, 6)
	226
	227
	228
	229	def test_boundary_allocation(self):
	230	quantity = Quantity(self.mesh4,
	231	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	232
	233	assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
	234
	235
	236	def test_set_values(self):
	237	quantity = Quantity(self.mesh4)
	238
	239
	240	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
	241	location = 'vertices')
	242	assert allclose(quantity.vertex_values,
	243	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	244	assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
	245	assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
	246	[5., 5., 5.],
	247	[4.5, 4.5, 0.],
	248	[3.0, -1.5, -1.5]])
	249
	250
	251	# Test default
	252	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	253	assert allclose(quantity.vertex_values,
	254	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	255	assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
	256	assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
	257	[5., 5., 5.],
	258	[4.5, 4.5, 0.],
	259	[3.0, -1.5, -1.5]])
	260
	261	# Test centroids
	262	quantity.set_values([1,2,3,4], location = 'centroids')
	263	assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
	264
	265	# Test exceptions
	266	try:
	267	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
	268	location = 'bas kamel tuba')
	269	except:
	270	pass
	271
	272
	273	try:
	274	quantity.set_values([[1,2,3], [0,0,9]])
	275	except AssertionError:
	276	pass
	277	except:
	278	raise 'should have raised Assertionerror'
	279
	280
	281
	282	def test_set_values_const(self):
	283	quantity = Quantity(self.mesh4)
	284
	285	quantity.set_values(1.0, location = 'vertices')
	286	assert allclose(quantity.vertex_values,
	287	[[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
	288
	289	assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
	290	assert allclose(quantity.edge_values, [[1, 1, 1],
	291	[1, 1, 1],
	292	[1, 1, 1],
	293	[1, 1, 1]])
	294
	295
	296	quantity.set_values(2.0, location = 'centroids')
	297	assert allclose(quantity.centroid_values, [2, 2, 2, 2])
	298
	299
	300	def test_set_values_func(self):
	301	quantity = Quantity(self.mesh4)
	302
	303	def f(x, y):
	304	return x+y
	305
	306	quantity.set_values(f, location = 'vertices')
	307	#print "quantity.vertex_values",quantity.vertex_values
	308	assert allclose(quantity.vertex_values,
	309	[[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
	310	assert allclose(quantity.centroid_values,
	311	[4.0/3, 8.0/3, 10.0/3, 10.0/3])
	312	assert allclose(quantity.edge_values,
	313	[[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
	314
	315
	316	quantity.set_values(f, location = 'centroids')
	317	assert allclose(quantity.centroid_values,
	318	[4.0/3, 8.0/3, 10.0/3, 10.0/3])
	319
	320
	321	def test_integral(self):
	322	quantity = Quantity(self.mesh4)
	323
	324	#Try constants first
	325	const = 5
	326	quantity.set_values(const, location = 'vertices')
	327	#print 'Q', quantity.get_integral()
	328
	329	assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)
	330
	331	#Try with a linear function
	332	def f(x, y):
	333	return x+y
	334
	335	quantity.set_values(f, location = 'vertices')
	336
	337
	338	ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
	339
	340	assert allclose (quantity.get_integral(), ref_integral)
	341
	342
	343
	344	def test_set_vertex_values(self):
	345	quantity = Quantity(self.mesh4)
	346	quantity.set_vertex_values([0,1,2,3,4,5])
	347
	348	assert allclose(quantity.vertex_values,
	349	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
	350	assert allclose(quantity.centroid_values,
	351	[1., 7./3, 11./3, 8./3]) #Centroid
	352	assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
	353	[3., 2.5, 1.5],
	354	[3.5, 4.5, 3.],
	355	[2.5, 3.5, 2]])
	356
	357
	358	def test_set_vertex_values_subset(self):
	359	quantity = Quantity(self.mesh4)
	360	quantity.set_vertex_values([0,1,2,3,4,5])
	361	quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
	362
	363	assert allclose(quantity.vertex_values,
	364	[[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
	365
	366
	367	def test_set_vertex_values_using_general_interface(self):
	368	quantity = Quantity(self.mesh4)
	369
	370
	371	quantity.set_values([0,1,2,3,4,5])
	372
	373
	374	assert allclose(quantity.vertex_values,
	375	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
	376
	377	#Centroid
	378	assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
	379
	380	assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
	381	[3., 2.5, 1.5],
	382	[3.5, 4.5, 3.],
	383	[2.5, 3.5, 2]])
	384
	385
	386
	387	def test_set_vertex_values_using_general_interface_with_subset(self):
	388	"""test_set_vertex_values_using_general_interface_with_subset(self):
	389
	390	Test that indices and polygon works (for constants values)
	391	"""
	392
	393	quantity = Quantity(self.mesh4)
	394
	395
	396	quantity.set_values([0,2,3,5], indices=[0,2,3,5])
	397	assert allclose(quantity.vertex_values,
	398	[[0,0,2], [0,2,0], [0,2,5], [3,0,0]])
	399
	400
	401	# Constant
	402	quantity.set_values(0.0)
	403	quantity.set_values(3.14, indices=[0,2], location='vertices')
	404
	405	# Indices refer to triangle numbers
	406	assert allclose(quantity.vertex_values,
	407	[[3.14,3.14,3.14], [0,0,0],
	408	[3.14,3.14,3.14], [0,0,0]])
	409
	410
	411
	412	# Now try with polygon (pick points where y>2)
	413	polygon = [[0,2.1], [4,2.1], [4,7], [0,7]]
	414	quantity.set_values(0.0)
	415	quantity.set_values(3.14, polygon=polygon)
	416
	417	assert allclose(quantity.vertex_values,
	418	[[0,0,0], [0,0,0], [0,0,0],
	419	[3.14,3.14,3.14]])
	420
	421
	422	# Another polygon (pick triangle 1 and 2 (rightmost triangles)
	423	# using centroids
	424	polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
	425	quantity.set_values(0.0)
	426	quantity.set_values(3.14, location='centroids', polygon=polygon)
	427	assert allclose(quantity.vertex_values,
	428	[[0,0,0],
	429	[3.14,3.14,3.14],
	430	[3.14,3.14,3.14],
	431	[0,0,0]])
	432
	433
	434	# Same polygon now use vertices (default)
	435	polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
	436	quantity.set_values(0.0)
	437	#print 'Here 2'
	438	quantity.set_values(3.14, polygon=polygon)
	439	assert allclose(quantity.vertex_values,
	440	[[0,0,0],
	441	[3.14,3.14,3.14],
	442	[3.14,3.14,3.14],
	443	[0,0,0]])
	444
	445
	446	# Test input checking
	447	try:
	448	quantity.set_values(3.14, polygon=polygon, indices = [0,2])
	449	except:
	450	pass
	451	else:
	452	msg = 'Should have caught this'
	453	raise msg
	454
	455
	456
	457
	458
	459	def test_set_vertex_values_using_general_interface_subset_and_geo(self):
	460	"""test_set_vertex_values_using_general_interface_with_subset(self):
	461	Test that indices and polygon works using georeferencing
	462	"""
	463
	464	quantity = Quantity(self.mesh4)
	465	G = Geo_reference(56, 10, 100)
	466	quantity.domain.geo_reference = G
	467
	468	#print quantity.domain.get_nodes(absolute=True)
	469
	470
	471	# Constant
	472	quantity.set_values(0.0)
	473	quantity.set_values(3.14, indices=[0,2], location='vertices')
	474
	475	# Indices refer to triangle numbers here - not vertices (why?)
	476	assert allclose(quantity.vertex_values,
	477	[[3.14,3.14,3.14], [0,0,0],
	478	[3.14,3.14,3.14], [0,0,0]])
	479
	480
	481
	482	# Now try with polygon (pick points where y>2)
	483	polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]])
	484	polygon += [G.xllcorner, G.yllcorner]
	485
	486	quantity.set_values(0.0)
	487	quantity.set_values(3.14, polygon=polygon, location='centroids')
	488
	489	assert allclose(quantity.vertex_values,
	490	[[0,0,0], [0,0,0], [0,0,0],
	491	[3.14,3.14,3.14]])
	492
	493
	494	# Another polygon (pick triangle 1 and 2 (rightmost triangles)
	495	polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]])
	496	polygon += [G.xllcorner, G.yllcorner]
	497
	498	quantity.set_values(0.0)
	499	quantity.set_values(3.14, polygon=polygon)
	500
	501	assert allclose(quantity.vertex_values,
	502	[[0,0,0],
	503	[3.14,3.14,3.14],
	504	[3.14,3.14,3.14],
	505	[0,0,0]])
	506
	507
	508
	509	def test_set_values_using_fit(self):
	510
	511
	512	quantity = Quantity(self.mesh4)
	513
	514	#Get (enough) datapoints
	515	data_points = [[ 0.66666667, 0.66666667],
	516	[ 1.33333333, 1.33333333],
	517	[ 2.66666667, 0.66666667],
	518	[ 0.66666667, 2.66666667],
	519	[ 0.0, 1.0],
	520	[ 0.0, 3.0],
	521	[ 1.0, 0.0],
	522	[ 1.0, 1.0],
	523	[ 1.0, 2.0],
	524	[ 1.0, 3.0],
	525	[ 2.0, 1.0],
	526	[ 3.0, 0.0],
	527	[ 3.0, 1.0]]
	528
	529	z = linear_function(data_points)
	530
	531	#Use built-in fit_interpolate.fit
	532	quantity.set_values( Geospatial_data(data_points, z), alpha = 0 )
	533	#quantity.set_values(points = data_points, values = z, alpha = 0)
	534
	535
	536	answer = linear_function(quantity.domain.get_vertex_coordinates())
	537	#print quantity.vertex_values, answer
	538	assert allclose(quantity.vertex_values.flat, answer)
	539
	540
	541	#Now try by setting the same values directly
	542	vertex_attributes = fit_to_mesh(data_points,
	543	quantity.domain.get_nodes(),
	544	quantity.domain.triangles, #FIXME
	545	point_attributes=z,
	546	alpha = 0,
	547	verbose=False)
	548
	549	#print vertex_attributes
	550	quantity.set_values(vertex_attributes)
	551	assert allclose(quantity.vertex_values.flat, answer)
	552
	553
	554
	555
	556
	557	def test_test_set_values_using_fit_w_geo(self):
	558
	559
	560	#Mesh
	561	vertex_coordinates = [[0.76, 0.76],
	562	[0.76, 5.76],
	563	[5.76, 0.76]]
	564	triangles = [[0,2,1]]
	565
	566	mesh_georef = Geo_reference(56,-0.76,-0.76)
	567	mesh1 = Domain(vertex_coordinates, triangles,
	568	geo_reference = mesh_georef)
	569	mesh1.check_integrity()
	570
	571	#Quantity
	572	quantity = Quantity(mesh1)
	573
	574	#Data
	575	data_points = [[ 201.0, 401.0],
	576	[ 201.0, 403.0],
	577	[ 203.0, 401.0]]
	578
	579	z = [2, 4, 4]
	580
	581	data_georef = Geo_reference(56,-200,-400)
	582
	583
	584	#Reference
	585	ref = fit_to_mesh(data_points, vertex_coordinates, triangles,
	586	point_attributes=z,
	587	data_origin = data_georef.get_origin(),
	588	mesh_origin = mesh_georef.get_origin(),
	589	alpha = 0)
	590
	591	assert allclose( ref, [0,5,5] )
	592
	593
	594	#Test set_values
	595
	596	quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 )
	597
	598	#quantity.set_values(points = data_points,
	599	# values = z,
	600	# data_georef = data_georef,
	601	# alpha = 0)
	602
	603
	604	#quantity.set_values(points = data_points,
	605	# values = z,
	606	# data_georef = data_georef,
	607	# alpha = 0)
	608	assert allclose(quantity.vertex_values.flat, ref)
	609
	610
	611
	612	#Test set_values using geospatial data object
	613	quantity.vertex_values[:] = 0.0
	614
	615	geo = Geospatial_data(data_points, z, data_georef)
	616
	617
	618	quantity.set_values(geospatial_data = geo, alpha = 0)
	619	assert allclose(quantity.vertex_values.flat, ref)
	620
	621
	622
	623	def test_set_values_from_file1(self):
	624	quantity = Quantity(self.mesh4)
	625
	626	#Get (enough) datapoints
	627	data_points = [[ 0.66666667, 0.66666667],
	628	[ 1.33333333, 1.33333333],
	629	[ 2.66666667, 0.66666667],
	630	[ 0.66666667, 2.66666667],
	631	[ 0.0, 1.0],
	632	[ 0.0, 3.0],
	633	[ 1.0, 0.0],
	634	[ 1.0, 1.0],
	635	[ 1.0, 2.0],
	636	[ 1.0, 3.0],
	637	[ 2.0, 1.0],
	638	[ 3.0, 0.0],
	639	[ 3.0, 1.0]]
	640
	641	data_geo_spatial = Geospatial_data(data_points,
	642	geo_reference = Geo_reference(56, 0, 0))
	643	data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
	644	attributes = linear_function(data_points_absolute)
	645	att = 'spam_and_eggs'
	646
	647	#Create .txt file
	648	ptsfile = tempfile.mktemp(".txt")
	649	file = open(ptsfile,"w")
	650	file.write(" x,y," + att + " \n")
	651	for data_point, attribute in map(None, data_points_absolute
	652	,attributes):
	653	row = str(data_point[0]) + ',' + str(data_point[1]) \
	654	+ ',' + str(attribute)
	655	file.write(row + "\n")
	656	file.close()
	657
	658
	659	#Check that values can be set from file
	660	quantity.set_values(filename = ptsfile,
	661	attribute_name = att, alpha = 0)
	662	answer = linear_function(quantity.domain.get_vertex_coordinates())
	663
	664	#print quantity.vertex_values.flat
	665	#print answer
	666
	667
	668	assert allclose(quantity.vertex_values.flat, answer)
	669
	670
	671	#Check that values can be set from file using default attribute
	672	quantity.set_values(filename = ptsfile, alpha = 0)
	673	assert allclose(quantity.vertex_values.flat, answer)
	674
	675	#Cleanup
	676	import os
	677	os.remove(ptsfile)
	678
	679
	680
	681	def Xtest_set_values_from_file_using_polygon(self):
	682	"""test_set_values_from_file_using_polygon(self):
	683
	684	Test that polygon restriction works for general points data
	685	"""
	686
	687	quantity = Quantity(self.mesh4)
	688
	689	#Get (enough) datapoints
	690	data_points = [[ 0.66666667, 0.66666667],
	691	[ 1.33333333, 1.33333333],
	692	[ 2.66666667, 0.66666667],
	693	[ 0.66666667, 2.66666667],
	694	[ 0.0, 1.0],
	695	[ 0.0, 3.0],
	696	[ 1.0, 0.0],
	697	[ 1.0, 1.0],
	698	[ 1.0, 2.0],
	699	[ 1.0, 3.0],
	700	[ 2.0, 1.0],
	701	[ 3.0, 0.0],
	702	[ 3.0, 1.0]]
	703
	704	data_geo_spatial = Geospatial_data(data_points,
	705	geo_reference = Geo_reference(56, 0, 0))
	706	data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
	707	attributes = linear_function(data_points_absolute)
	708	att = 'spam_and_eggs'
	709
	710	#Create .txt file
	711	ptsfile = tempfile.mktemp(".txt")
	712	file = open(ptsfile,"w")
	713	file.write(" x,y," + att + " \n")
	714	for data_point, attribute in map(None, data_points_absolute
	715	,attributes):
	716	row = str(data_point[0]) + ',' + str(data_point[1]) \
	717	+ ',' + str(attribute)
	718	file.write(row + "\n")
	719	file.close()
	720
	721	# Create restricting polygon (containing node #4 (2,2) and
	722	# centroid of triangle #1 (bce)
	723	polygon = [[1.0, 1.0], [4.0, 1.0],
	724	[4.0, 4.0], [1.0, 4.0]]
	725
	726	#print self.mesh4.nodes
	727	#print inside_polygon(self.mesh4.nodes, polygon)
	728	assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)
	729
	730	#print quantity.domain.get_vertex_coordinates()
	731	#print quantity.domain.get_nodes()
	732
	733	# Check that values can be set from file
	734	quantity.set_values(filename=ptsfile,
	735	polygon=polygon,
	736	location='unique vertices',
	737	alpha=0)
	738
	739	# Get indices for vertex coordinates in polygon
	740	indices = inside_polygon(quantity.domain.get_vertex_coordinates(),
	741	polygon)
	742	points = take(quantity.domain.get_vertex_coordinates(), indices)
	743
	744	answer = linear_function(points)
	745
	746	#print quantity.vertex_values.flat
	747	#print answer
	748
	749	# Check vertices in polygon have been set
	750	assert allclose(take(quantity.vertex_values.flat, indices),
	751	answer)
	752
	753	# Check vertices outside polygon are zero
	754	indices = outside_polygon(quantity.domain.get_vertex_coordinates(),
	755	polygon)
	756	assert allclose(take(quantity.vertex_values.flat, indices),
	757	0.0)
	758
	759	#Cleanup
	760	import os
	761	os.remove(ptsfile)
	762
	763
	764
	765
	766	def test_cache_test_set_values_from_file(self):
	767	# FIXME (Ole): What is this about?
	768	# I don't think it checks anything new
	769	quantity = Quantity(self.mesh4)
	770
	771	#Get (enough) datapoints
	772	data_points = [[ 0.66666667, 0.66666667],
	773	[ 1.33333333, 1.33333333],
	774	[ 2.66666667, 0.66666667],
	775	[ 0.66666667, 2.66666667],
	776	[ 0.0, 1.0],
	777	[ 0.0, 3.0],
	778	[ 1.0, 0.0],
	779	[ 1.0, 1.0],
	780	[ 1.0, 2.0],
	781	[ 1.0, 3.0],
	782	[ 2.0, 1.0],
	783	[ 3.0, 0.0],
	784	[ 3.0, 1.0]]
	785
	786	georef = Geo_reference(56, 0, 0)
	787	data_geo_spatial = Geospatial_data(data_points,
	788	geo_reference=georef)
	789
	790	data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
	791	attributes = linear_function(data_points_absolute)
	792	att = 'spam_and_eggs'
	793
	794	# Create .txt file
	795	ptsfile = tempfile.mktemp(".txt")
	796	file = open(ptsfile,"w")
	797	file.write(" x,y," + att + " \n")
	798	for data_point, attribute in map(None, data_points_absolute
	799	,attributes):
	800	row = str(data_point[0]) + ',' + str(data_point[1]) \
	801	+ ',' + str(attribute)
	802	file.write(row + "\n")
	803	file.close()
	804
	805
	806	# Check that values can be set from file
	807	quantity.set_values(filename=ptsfile,
	808	attribute_name=att,
	809	alpha=0,
	810	use_cache=True,
	811	verbose=False)
	812	answer = linear_function(quantity.domain.get_vertex_coordinates())
	813	assert allclose(quantity.vertex_values.flat, answer)
	814
	815
	816	# Check that values can be set from file using default attribute
	817	quantity.set_values(filename=ptsfile,
	818	alpha=0)
	819	assert allclose(quantity.vertex_values.flat, answer)
	820
	821	# Check cache
	822	quantity.set_values(filename=ptsfile,
	823	attribute_name=att,
	824	alpha=0,
	825	use_cache=True,
	826	verbose=False)
	827
	828
	829	#Cleanup
	830	import os
	831	os.remove(ptsfile)
	832
	833	def test_set_values_from_lat_long(self):
	834	quantity = Quantity(self.mesh_onslow)
	835
	836	#Get (enough) datapoints
	837	data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	838	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
	839
	840	data_geo_spatial = Geospatial_data(data_points,
	841	points_are_lats_longs=True)
	842	points_UTM = data_geo_spatial.get_data_points(absolute=True)
	843	attributes = linear_function(points_UTM)
	844	att = 'elevation'
	845
	846	#Create .txt file
	847	txt_file = tempfile.mktemp(".txt")
	848	file = open(txt_file,"w")
	849	file.write(" lat,long," + att + " \n")
	850	for data_point, attribute in map(None, data_points, attributes):
	851	row = str(data_point[0]) + ',' + str(data_point[1]) \
	852	+ ',' + str(attribute)
	853	#print "row", row
	854	file.write(row + "\n")
	855	file.close()
	856
	857
	858	#Check that values can be set from file
	859	quantity.set_values(filename=txt_file,
	860	attribute_name=att,
	861	alpha=0)
	862	answer = linear_function(quantity.domain.get_vertex_coordinates())
	863
	864	#print "quantity.vertex_values.flat", quantity.vertex_values.flat
	865	#print "answer",answer
	866
	867	assert allclose(quantity.vertex_values.flat, answer)
	868
	869
	870	#Check that values can be set from file using default attribute
	871	quantity.set_values(filename=txt_file, alpha=0)
	872	assert allclose(quantity.vertex_values.flat, answer)
	873
	874	#Cleanup
	875	import os
	876	os.remove(txt_file)
	877
	878	def test_set_values_from_lat_long(self):
	879	quantity = Quantity(self.mesh_onslow)
	880
	881	#Get (enough) datapoints
	882	data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	883	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
	884
	885	data_geo_spatial = Geospatial_data(data_points,
	886	points_are_lats_longs=True)
	887	points_UTM = data_geo_spatial.get_data_points(absolute=True)
	888	attributes = linear_function(points_UTM)
	889	att = 'elevation'
	890
	891	#Create .txt file
	892	txt_file = tempfile.mktemp(".txt")
	893	file = open(txt_file,"w")
	894	file.write(" lat,long," + att + " \n")
	895	for data_point, attribute in map(None, data_points, attributes):
	896	row = str(data_point[0]) + ',' + str(data_point[1]) \
	897	+ ',' + str(attribute)
	898	#print "row", row
	899	file.write(row + "\n")
	900	file.close()
	901
	902
	903	#Check that values can be set from file
	904	quantity.set_values(filename=txt_file,
	905	attribute_name=att, alpha=0)
	906	answer = linear_function(quantity.domain.get_vertex_coordinates())
	907
	908	#print "quantity.vertex_values.flat", quantity.vertex_values.flat
	909	#print "answer",answer
	910
	911	assert allclose(quantity.vertex_values.flat, answer)
	912
	913
	914	#Check that values can be set from file using default attribute
	915	quantity.set_values(filename=txt_file, alpha=0)
	916	assert allclose(quantity.vertex_values.flat, answer)
	917
	918	#Cleanup
	919	import os
	920	os.remove(txt_file)
	921
	922	def test_set_values_from_UTM_pts(self):
	923	quantity = Quantity(self.mesh_onslow)
	924
	925	#Get (enough) datapoints
	926	data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	927	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
	928
	929	data_geo_spatial = Geospatial_data(data_points,
	930	points_are_lats_longs=True)
	931	points_UTM = data_geo_spatial.get_data_points(absolute=True)
	932	attributes = linear_function(points_UTM)
	933	att = 'elevation'
	934
	935	#Create .txt file
	936	txt_file = tempfile.mktemp(".txt")
	937	file = open(txt_file,"w")
	938	file.write(" x,y," + att + " \n")
	939	for data_point, attribute in map(None, points_UTM, attributes):
	940	row = str(data_point[0]) + ',' + str(data_point[1]) \
	941	+ ',' + str(attribute)
	942	#print "row", row
	943	file.write(row + "\n")
	944	file.close()
	945
	946
	947	pts_file = tempfile.mktemp(".pts")
	948	convert = Geospatial_data(txt_file)
	949	convert.export_points_file(pts_file)
	950
	951	#Check that values can be set from file
	952	quantity.set_values_from_file(pts_file, att, 0,
	953	'vertices', None)
	954	answer = linear_function(quantity.domain.get_vertex_coordinates())
	955	#print "quantity.vertex_values.flat", quantity.vertex_values.flat
	956	#print "answer",answer
	957	assert allclose(quantity.vertex_values.flat, answer)
	958
	959	#Check that values can be set from file
	960	quantity.set_values(filename=pts_file,
	961	attribute_name=att, alpha=0)
	962	answer = linear_function(quantity.domain.get_vertex_coordinates())
	963	#print "quantity.vertex_values.flat", quantity.vertex_values.flat
	964	#print "answer",answer
	965	assert allclose(quantity.vertex_values.flat, answer)
	966
	967
	968	#Check that values can be set from file using default attribute
	969	quantity.set_values(filename=txt_file, alpha=0)
	970	assert allclose(quantity.vertex_values.flat, answer)
	971
	972	#Cleanup
	973	import os
	974	os.remove(txt_file)
	975	os.remove(pts_file)
	976
	977	def verbose_test_set_values_from_UTM_pts(self):
	978	quantity = Quantity(self.mesh_onslow)
	979
	980	#Get (enough) datapoints
	981	data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	982	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	983	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	984	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	985	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	986	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	987	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	988	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	989	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	990	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	991	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	992	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	993	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	994	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	995	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	996	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	997	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	998	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	999	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	1000	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	1001	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	1002	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	1003	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	1004	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	1005	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	1006	[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
	1007	[-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
	1008	]
	1009
	1010	data_geo_spatial = Geospatial_data(data_points,
	1011	points_are_lats_longs=True)
	1012	points_UTM = data_geo_spatial.get_data_points(absolute=True)
	1013	attributes = linear_function(points_UTM)
	1014	att = 'elevation'
	1015
	1016	#Create .txt file
	1017	txt_file = tempfile.mktemp(".txt")
	1018	file = open(txt_file,"w")
	1019	file.write(" x,y," + att + " \n")
	1020	for data_point, attribute in map(None, points_UTM, attributes):
	1021	row = str(data_point[0]) + ',' + str(data_point[1]) \
	1022	+ ',' + str(attribute)
	1023	#print "row", row
	1024	file.write(row + "\n")
	1025	file.close()
	1026
	1027
	1028	pts_file = tempfile.mktemp(".pts")
	1029	convert = Geospatial_data(txt_file)
	1030	convert.export_points_file(pts_file)
	1031
	1032	#Check that values can be set from file
	1033	quantity.set_values_from_file(pts_file, att, 0,
	1034	'vertices', None, verbose = True,
	1035	max_read_lines=2)
	1036	answer = linear_function(quantity.domain.get_vertex_coordinates())
	1037	#print "quantity.vertex_values.flat", quantity.vertex_values.flat
	1038	#print "answer",answer
	1039	assert allclose(quantity.vertex_values.flat, answer)
	1040
	1041	#Check that values can be set from file
	1042	quantity.set_values(filename=pts_file,
	1043	attribute_name=att, alpha=0)
	1044	answer = linear_function(quantity.domain.get_vertex_coordinates())
	1045	#print "quantity.vertex_values.flat", quantity.vertex_values.flat
	1046	#print "answer",answer
	1047	assert allclose(quantity.vertex_values.flat, answer)
	1048
	1049
	1050	#Check that values can be set from file using default attribute
	1051	quantity.set_values(filename=txt_file, alpha=0)
	1052	assert allclose(quantity.vertex_values.flat, answer)
	1053
	1054	#Cleanup
	1055	import os
	1056	os.remove(txt_file)
	1057	os.remove(pts_file)
	1058
	1059	def test_set_values_from_file_with_georef1(self):
	1060
	1061	#Mesh in zone 56 (absolute coords)
	1062
	1063	x0 = 314036.58727982
	1064	y0 = 6224951.2960092
	1065
	1066	a = [x0+0.0, y0+0.0]
	1067	b = [x0+0.0, y0+2.0]
	1068	c = [x0+2.0, y0+0.0]
	1069	d = [x0+0.0, y0+4.0]
	1070	e = [x0+2.0, y0+2.0]
	1071	f = [x0+4.0, y0+0.0]
	1072
	1073	points = [a, b, c, d, e, f]
	1074
	1075	#bac, bce, ecf, dbe
	1076	elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
	1077
	1078	#absolute going in ..
	1079	mesh4 = Domain(points, elements,
	1080	geo_reference = Geo_reference(56, 0, 0))
	1081	mesh4.check_integrity()
	1082	quantity = Quantity(mesh4)
	1083
	1084	#Get (enough) datapoints (relative to georef)
	1085	data_points_rel = [[ 0.66666667, 0.66666667],
	1086	[ 1.33333333, 1.33333333],
	1087	[ 2.66666667, 0.66666667],
	1088	[ 0.66666667, 2.66666667],
	1089	[ 0.0, 1.0],
	1090	[ 0.0, 3.0],
	1091	[ 1.0, 0.0],
	1092	[ 1.0, 1.0],
	1093	[ 1.0, 2.0],
	1094	[ 1.0, 3.0],
	1095	[ 2.0, 1.0],
	1096	[ 3.0, 0.0],
	1097	[ 3.0, 1.0]]
	1098
	1099	data_geo_spatial = Geospatial_data(data_points_rel,
	1100	geo_reference = Geo_reference(56, x0, y0))
	1101	data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
	1102	attributes = linear_function(data_points_absolute)
	1103	att = 'spam_and_eggs'
	1104
	1105	#Create .txt file
	1106	ptsfile = tempfile.mktemp(".txt")
	1107	file = open(ptsfile,"w")
	1108	file.write(" x,y," + att + " \n")
	1109	for data_point, attribute in map(None, data_points_absolute
	1110	,attributes):
	1111	row = str(data_point[0]) + ',' + str(data_point[1]) \
	1112	+ ',' + str(attribute)
	1113	file.write(row + "\n")
	1114	file.close()
	1115
	1116	#file = open(ptsfile, 'r')
	1117	#lines = file.readlines()
	1118	#file.close()
	1119
	1120
	1121	#Check that values can be set from file
	1122	quantity.set_values(filename=ptsfile,
	1123	attribute_name=att, alpha=0)
	1124	answer = linear_function(quantity.domain.get_vertex_coordinates())
	1125
	1126	assert allclose(quantity.vertex_values.flat, answer)
	1127
	1128
	1129	#Check that values can be set from file using default attribute
	1130	quantity.set_values(filename=ptsfile, alpha=0)
	1131	assert allclose(quantity.vertex_values.flat, answer)
	1132
	1133	#Cleanup
	1134	import os
	1135	os.remove(ptsfile)
	1136
	1137
	1138	def test_set_values_from_file_with_georef2(self):
	1139
	1140	#Mesh in zone 56 (relative coords)
	1141
	1142	x0 = 314036.58727982
	1143	y0 = 6224951.2960092
	1144	#x0 = 0.0
	1145	#y0 = 0.0
	1146
	1147	a = [0.0, 0.0]
	1148	b = [0.0, 2.0]
	1149	c = [2.0, 0.0]
	1150	d = [0.0, 4.0]
	1151	e = [2.0, 2.0]
	1152	f = [4.0, 0.0]
	1153
	1154	points = [a, b, c, d, e, f]
	1155
	1156	#bac, bce, ecf, dbe
	1157	elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
	1158
	1159	mesh4 = Domain(points, elements,
	1160	geo_reference = Geo_reference(56, x0, y0))
	1161	mesh4.check_integrity()
	1162	quantity = Quantity(mesh4)
	1163
	1164	#Get (enough) datapoints
	1165	data_points = [[ x0+0.66666667, y0+0.66666667],
	1166	[ x0+1.33333333, y0+1.33333333],
	1167	[ x0+2.66666667, y0+0.66666667],
	1168	[ x0+0.66666667, y0+2.66666667],
	1169	[ x0+0.0, y0+1.0],
	1170	[ x0+0.0, y0+3.0],
	1171	[ x0+1.0, y0+0.0],
	1172	[ x0+1.0, y0+1.0],
	1173	[ x0+1.0, y0+2.0],
	1174	[ x0+1.0, y0+3.0],
	1175	[ x0+2.0, y0+1.0],
	1176	[ x0+3.0, y0+0.0],
	1177	[ x0+3.0, y0+1.0]]
	1178
	1179
	1180	data_geo_spatial = Geospatial_data(data_points,
	1181	geo_reference = Geo_reference(56, 0, 0))
	1182	data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
	1183	attributes = linear_function(data_points_absolute)
	1184	att = 'spam_and_eggs'
	1185
	1186	#Create .txt file
	1187	ptsfile = tempfile.mktemp(".txt")
	1188	file = open(ptsfile,"w")
	1189	file.write(" x,y," + att + " \n")
	1190	for data_point, attribute in map(None, data_points_absolute
	1191	,attributes):
	1192	row = str(data_point[0]) + ',' + str(data_point[1]) \
	1193	+ ',' + str(attribute)
	1194	file.write(row + "\n")
	1195	file.close()
	1196
	1197
	1198	#Check that values can be set from file
	1199	quantity.set_values(filename=ptsfile,
	1200	attribute_name=att, alpha=0)
	1201	answer = linear_function(quantity.domain. \
	1202	get_vertex_coordinates(absolute=True))
	1203
	1204
	1205	assert allclose(quantity.vertex_values.flat, answer)
	1206
	1207
	1208	#Check that values can be set from file using default attribute
	1209	quantity.set_values(filename=ptsfile, alpha=0)
	1210	assert allclose(quantity.vertex_values.flat, answer)
	1211
	1212	#Cleanup
	1213	import os
	1214	os.remove(ptsfile)
	1215
	1216
	1217
	1218
	1219	def test_set_values_from_quantity(self):
	1220
	1221	quantity1 = Quantity(self.mesh4)
	1222	quantity1.set_vertex_values([0,1,2,3,4,5])
	1223
	1224	assert allclose(quantity1.vertex_values,
	1225	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
	1226
	1227
	1228	quantity2 = Quantity(self.mesh4)
	1229	quantity2.set_values(quantity=quantity1)
	1230	assert allclose(quantity2.vertex_values,
	1231	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
	1232
	1233	quantity2.set_values(quantity = 2*quantity1)
	1234	assert allclose(quantity2.vertex_values,
	1235	[[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
	1236
	1237	quantity2.set_values(quantity = 2*quantity1 + 3)
	1238	assert allclose(quantity2.vertex_values,
	1239	[[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
	1240
	1241
	1242	#Check detection of quantity as first orgument
	1243	quantity2.set_values(2*quantity1 + 3)
	1244	assert allclose(quantity2.vertex_values,
	1245	[[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
	1246
	1247
	1248
	1249	def Xtest_set_values_from_quantity_using_polygon(self):
	1250	"""test_set_values_from_quantity_using_polygon(self):
	1251
	1252	Check that polygon can be used to restrict set_values when
	1253	using another quantity as argument.
	1254	"""
	1255
	1256	# Create restricting polygon (containing node #4 (2,2) and
	1257	# centroid of triangle #1 (bce)
	1258	polygon = [[1.0, 1.0], [4.0, 1.0],
	1259	[4.0, 4.0], [1.0, 4.0]]
	1260	assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)
	1261
	1262	quantity1 = Quantity(self.mesh4)
	1263	quantity1.set_vertex_values([0,1,2,3,4,5])
	1264
	1265	assert allclose(quantity1.vertex_values,
	1266	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
	1267
	1268
	1269	quantity2 = Quantity(self.mesh4)
	1270	quantity2.set_values(quantity=quantity1,
	1271	polygon=polygon)
	1272
	1273	msg = 'Only node #4(e) at (2,2) should have values applied '
	1274	assert allclose(quantity2.vertex_values,
	1275	[[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg
	1276	#bac, bce, ecf, dbe
	1277
	1278
	1279
	1280	def test_overloading(self):
	1281
	1282	quantity1 = Quantity(self.mesh4)
	1283	quantity1.set_vertex_values([0,1,2,3,4,5])
	1284
	1285	assert allclose(quantity1.vertex_values,
	1286	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
	1287
	1288
	1289	quantity2 = Quantity(self.mesh4)
	1290	quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
	1291	location = 'vertices')
	1292
	1293
	1294
	1295	quantity3 = Quantity(self.mesh4)
	1296	quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]],
	1297	location = 'vertices')
	1298
	1299
	1300	# Negation
	1301	Q = -quantity1
	1302	assert allclose(Q.vertex_values, -quantity1.vertex_values)
	1303	assert allclose(Q.centroid_values, -quantity1.centroid_values)
	1304	assert allclose(Q.edge_values, -quantity1.edge_values)
	1305
	1306	# Addition
	1307	Q = quantity1 + 7
	1308	assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
	1309	assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
	1310	assert allclose(Q.edge_values, quantity1.edge_values + 7)
	1311
	1312	Q = 7 + quantity1
	1313	assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
	1314	assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
	1315	assert allclose(Q.edge_values, quantity1.edge_values + 7)
	1316
	1317	Q = quantity1 + quantity2
	1318	assert allclose(Q.vertex_values,
	1319	quantity1.vertex_values + quantity2.vertex_values)
	1320	assert allclose(Q.centroid_values,
	1321	quantity1.centroid_values + quantity2.centroid_values)
	1322	assert allclose(Q.edge_values,
	1323	quantity1.edge_values + quantity2.edge_values)
	1324
	1325
	1326	Q = quantity1 + quantity2 - 3
	1327	assert allclose(Q.vertex_values,
	1328	quantity1.vertex_values + quantity2.vertex_values - 3)
	1329
	1330	Q = quantity1 - quantity2
	1331	assert allclose(Q.vertex_values,
	1332	quantity1.vertex_values - quantity2.vertex_values)
	1333
	1334	#Scaling
	1335	Q = quantity1*3
	1336	assert allclose(Q.vertex_values, quantity1.vertex_values*3)
	1337	assert allclose(Q.centroid_values, quantity1.centroid_values*3)
	1338	assert allclose(Q.edge_values, quantity1.edge_values*3)
	1339	Q = 3*quantity1
	1340	assert allclose(Q.vertex_values, quantity1.vertex_values*3)
	1341
	1342	#Multiplication
	1343	Q = quantity1 * quantity2
	1344	#print Q.vertex_values
	1345	#print Q.centroid_values
	1346	#print quantity1.centroid_values
	1347	#print quantity2.centroid_values
	1348
	1349	assert allclose(Q.vertex_values,
	1350	quantity1.vertex_values * quantity2.vertex_values)
	1351
	1352	#Linear combinations
	1353	Q = 4*quantity1 + 2
	1354	assert allclose(Q.vertex_values,
	1355	4*quantity1.vertex_values + 2)
	1356
	1357	Q = quantity1*quantity2 + 2
	1358	assert allclose(Q.vertex_values,
	1359	quantity1.vertex_values * quantity2.vertex_values + 2)
	1360
	1361	Q = quantity1*quantity2 + quantity3
	1362	assert allclose(Q.vertex_values,
	1363	quantity1.vertex_values * quantity2.vertex_values +
	1364	quantity3.vertex_values)
	1365	Q = quantity1quantity2 + 3quantity3
	1366	assert allclose(Q.vertex_values,
	1367	quantity1.vertex_values * quantity2.vertex_values +
	1368	3*quantity3.vertex_values)
	1369	Q = quantity1quantity2 + 3quantity3 + 5.0
	1370	assert allclose(Q.vertex_values,
	1371	quantity1.vertex_values * quantity2.vertex_values +
	1372	3*quantity3.vertex_values + 5)
	1373
	1374	Q = quantity1*quantity2 - quantity3
	1375	assert allclose(Q.vertex_values,
	1376	quantity1.vertex_values * quantity2.vertex_values -
	1377	quantity3.vertex_values)
	1378	Q = 1.5quantity1quantity2 - 3*quantity3 + 5.0
	1379	assert allclose(Q.vertex_values,
	1380	1.5quantity1.vertex_values quantity2.vertex_values -
	1381	3*quantity3.vertex_values + 5)
	1382
	1383	#Try combining quantities and arrays and scalars
	1384	Q = 1.5quantity1quantity2.vertex_values -\
	1385	3*quantity3.vertex_values + 5.0
	1386	assert allclose(Q.vertex_values,
	1387	1.5quantity1.vertex_values quantity2.vertex_values -
	1388	3*quantity3.vertex_values + 5)
	1389
	1390
	1391	#Powers
	1392	Q = quantity1**2
	1393	assert allclose(Q.vertex_values, quantity1.vertex_values**2)
	1394
	1395	Q = quantity12 +quantity22
	1396	assert allclose(Q.vertex_values,
	1397	quantity1.vertex_values**2 + \
	1398	quantity2.vertex_values**2)
	1399
	1400	Q = (quantity12 +quantity22)**0.5
	1401	assert allclose(Q.vertex_values,
	1402	(quantity1.vertex_values**2 + \
	1403	quantity2.vertex_values2)0.5)
	1404
	1405
	1406
	1407
	1408
	1409
	1410
	1411	def test_compute_gradient(self):
	1412	quantity = Quantity(self.mesh4)
	1413
	1414	#Set up for a gradient of (2,0) at mid triangle
	1415	quantity.set_values([2.0, 4.0, 6.0, 2.0],
	1416	location = 'centroids')
	1417
	1418
	1419	#Gradients
	1420	quantity.compute_gradients()
	1421
	1422	a = quantity.x_gradient
	1423	b = quantity.y_gradient
	1424	#print self.mesh4.centroid_coordinates
	1425	#print a, b
	1426
	1427	#The central triangle (1)
	1428	#(using standard gradient based on neigbours controid values)
	1429	assert allclose(a[1], 2.0)
	1430	assert allclose(b[1], 0.0)
	1431
	1432
	1433	#Left triangle (0) using two point gradient
	1434	#q0 = q1 + a(x0-x1) + b(y0-y1) <=>
	1435	#2 = 4 + a(-2/3) + b(-2/3)
	1436	assert allclose(a[0] + b[0], 3)
	1437	#From orthogonality (a(y0-y1) + b(x0-x1) == 0)
	1438	assert allclose(a[0] - b[0], 0)
	1439
	1440
	1441	#Right triangle (2) using two point gradient
	1442	#q2 = q1 + a(x2-x1) + b(y2-y1) <=>
	1443	#6 = 4 + a(4/3) + b(-2/3)
	1444	assert allclose(2*a[2] - b[2], 3)
	1445	#From orthogonality (a(y1-y2) + b(x2-x1) == 0)
	1446	assert allclose(a[2] + 2*b[2], 0)
	1447
	1448
	1449	#Top triangle (3) using two point gradient
	1450	#q3 = q1 + a(x3-x1) + b(y3-y1) <=>
	1451	#2 = 4 + a(-2/3) + b(4/3)
	1452	assert allclose(a[3] - 2*b[3], 3)
	1453	#From orthogonality (a(y1-y3) + b(x3-x1) == 0)
	1454	assert allclose(2*a[3] + b[3], 0)
	1455
	1456
	1457
	1458	#print a, b
	1459	quantity.extrapolate_second_order()
	1460
	1461	#Apply q(x,y) = qc + a(x-xc) + b(y-yc)
	1462	assert allclose(quantity.vertex_values[0,:], [3., 0., 3.])
	1463	assert allclose(quantity.vertex_values[1,:], [4./3, 16./3, 16./3])
	1464
	1465
	1466	#a = 1.2, b=-0.6
	1467	#q(4,0) = 6 + a(4 - 8/3) + b(-2/3)
	1468	assert allclose(quantity.vertex_values[2,2], 8)
	1469
	1470	def test_get_gradients(self):
	1471	quantity = Quantity(self.mesh4)
	1472
	1473	#Set up for a gradient of (2,0) at mid triangle
	1474	quantity.set_values([2.0, 4.0, 6.0, 2.0],
	1475	location = 'centroids')
	1476
	1477
	1478	#Gradients
	1479	quantity.compute_gradients()
	1480
	1481	a, b = quantity.get_gradients()
	1482	#print self.mesh4.centroid_coordinates
	1483	#print a, b
	1484
	1485	#The central triangle (1)
	1486	#(using standard gradient based on neigbours controid values)
	1487	assert allclose(a[1], 2.0)
	1488	assert allclose(b[1], 0.0)
	1489
	1490
	1491	#Left triangle (0) using two point gradient
	1492	#q0 = q1 + a(x0-x1) + b(y0-y1) <=>
	1493	#2 = 4 + a(-2/3) + b(-2/3)
	1494	assert allclose(a[0] + b[0], 3)
	1495	#From orthogonality (a(y0-y1) + b(x0-x1) == 0)
	1496	assert allclose(a[0] - b[0], 0)
	1497
	1498
	1499	#Right triangle (2) using two point gradient
	1500	#q2 = q1 + a(x2-x1) + b(y2-y1) <=>
	1501	#6 = 4 + a(4/3) + b(-2/3)
	1502	assert allclose(2*a[2] - b[2], 3)
	1503	#From orthogonality (a(y1-y2) + b(x2-x1) == 0)
	1504	assert allclose(a[2] + 2*b[2], 0)
	1505
	1506
	1507	#Top triangle (3) using two point gradient
	1508	#q3 = q1 + a(x3-x1) + b(y3-y1) <=>
	1509	#2 = 4 + a(-2/3) + b(4/3)
	1510	assert allclose(a[3] - 2*b[3], 3)
	1511	#From orthogonality (a(y1-y3) + b(x3-x1) == 0)
	1512	assert allclose(2*a[3] + b[3], 0)
	1513
	1514
	1515	def test_second_order_extrapolation2(self):
	1516	quantity = Quantity(self.mesh4)
	1517
	1518	#Set up for a gradient of (3,1), f(x) = 3x+y
	1519	quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
	1520	location = 'centroids')
	1521
	1522	#Gradients
	1523	quantity.compute_gradients()
	1524
	1525	a = quantity.x_gradient
	1526	b = quantity.y_gradient
	1527
	1528	#print a, b
	1529
	1530	assert allclose(a[1], 3.0)
	1531	assert allclose(b[1], 1.0)
	1532
	1533	#Work out the others
	1534
	1535	quantity.extrapolate_second_order()
	1536
	1537	#print quantity.vertex_values
	1538	assert allclose(quantity.vertex_values[1,0], 2.0)
	1539	assert allclose(quantity.vertex_values[1,1], 6.0)
	1540	assert allclose(quantity.vertex_values[1,2], 8.0)
	1541
	1542
	1543
	1544	def test_backup_saxpy_centroid_values(self):
	1545	quantity = Quantity(self.mesh4)
	1546
	1547	#Set up for a gradient of (3,1), f(x) = 3x+y
	1548	c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
	1549	d_values = array([1.0, 2.0, 3.0, 4.0])
	1550	quantity.set_values(c_values, location = 'centroids')
	1551
	1552	#Backup
	1553	quantity.backup_centroid_values()
	1554
	1555	#print quantity.vertex_values
	1556	assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
	1557
	1558
	1559	quantity.set_values(d_values, location = 'centroids')
	1560
	1561	quantity.saxpy_centroid_values(2.0, 3.0)
	1562
	1563	assert(quantity.centroid_values, 2.0d_values + 3.0c_values)
	1564
	1565
	1566
	1567	def test_first_order_extrapolator(self):
	1568	quantity = Quantity(self.mesh4)
	1569
	1570	#Test centroids
	1571	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
	1572	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
	1573
	1574	#Extrapolate
	1575	quantity.extrapolate_first_order()
	1576
	1577	#Check that gradient is zero
	1578	a,b = quantity.get_gradients()
	1579	assert allclose(a, [0,0,0,0])
	1580	assert allclose(b, [0,0,0,0])
	1581
	1582	#Check vertices but not edge values
	1583	assert allclose(quantity.vertex_values,
	1584	[[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
	1585
	1586
	1587	def test_second_order_extrapolator(self):
	1588	quantity = Quantity(self.mesh4)
	1589
	1590	#Set up for a gradient of (3,0) at mid triangle
	1591	quantity.set_values([2.0, 4.0, 8.0, 2.0],
	1592	location = 'centroids')
	1593
	1594
	1595
	1596	quantity.extrapolate_second_order()
	1597	quantity.limit()
	1598
	1599
	1600	#Assert that central triangle is limited by neighbours
	1601	assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
	1602	assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
	1603
	1604	assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
	1605	assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
	1606
	1607	assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
	1608	assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
	1609
	1610
	1611	#Assert that quantities are conserved
	1612	from Numeric import sum
	1613	for k in range(quantity.centroid_values.shape[0]):
	1614	assert allclose (quantity.centroid_values[k],
	1615	sum(quantity.vertex_values[k,:])/3)
	1616
	1617
	1618
	1619
	1620
	1621	def test_limit_vertices_by_all_neighbours(self):
	1622	quantity = Quantity(self.mesh4)
	1623
	1624	#Create a deliberate overshoot (e.g. from gradient computation)
	1625	quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
	1626
	1627
	1628	#Limit
	1629	quantity.limit_vertices_by_all_neighbours()
	1630
	1631	#Assert that central triangle is limited by neighbours
	1632	assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
	1633	assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
	1634
	1635	assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
	1636	assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
	1637
	1638	assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
	1639	assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
	1640
	1641
	1642
	1643	#Assert that quantities are conserved
	1644	from Numeric import sum
	1645	for k in range(quantity.centroid_values.shape[0]):
	1646	assert allclose (quantity.centroid_values[k],
	1647	sum(quantity.vertex_values[k,:])/3)
	1648
	1649
	1650
	1651	def test_limit_edges_by_all_neighbours(self):
	1652	quantity = Quantity(self.mesh4)
	1653
	1654	#Create a deliberate overshoot (e.g. from gradient computation)
	1655	quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
	1656
	1657
	1658	#Limit
	1659	quantity.limit_edges_by_all_neighbours()
	1660
	1661	#Assert that central triangle is limited by neighbours
	1662	assert quantity.edge_values[1,0] <= quantity.centroid_values[2]
	1663	assert quantity.edge_values[1,0] >= quantity.centroid_values[0]
	1664
	1665	assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
	1666	assert quantity.edge_values[1,1] >= quantity.centroid_values[0]
	1667
	1668	assert quantity.edge_values[1,2] <= quantity.centroid_values[2]
	1669	assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
	1670
	1671
	1672
	1673	#Assert that quantities are conserved
	1674	from Numeric import sum
	1675	for k in range(quantity.centroid_values.shape[0]):
	1676	assert allclose (quantity.centroid_values[k],
	1677	sum(quantity.vertex_values[k,:])/3)
	1678
	1679
	1680	def test_limit_edges_by_neighbour(self):
	1681	quantity = Quantity(self.mesh4)
	1682
	1683	#Create a deliberate overshoot (e.g. from gradient computation)
	1684	quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
	1685
	1686
	1687	#Limit
	1688	quantity.limit_edges_by_neighbour()
	1689
	1690	#Assert that central triangle is limited by neighbours
	1691	assert quantity.edge_values[1,0] <= quantity.centroid_values[3]
	1692	assert quantity.edge_values[1,0] >= quantity.centroid_values[1]
	1693
	1694	assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
	1695	assert quantity.edge_values[1,1] >= quantity.centroid_values[1]
	1696
	1697	assert quantity.edge_values[1,2] <= quantity.centroid_values[1]
	1698	assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
	1699
	1700
	1701
	1702	#Assert that quantities are conserved
	1703	from Numeric import sum
	1704	for k in range(quantity.centroid_values.shape[0]):
	1705	assert allclose (quantity.centroid_values[k],
	1706	sum(quantity.vertex_values[k,:])/3)
	1707
	1708	def test_limiter2(self):
	1709	"""Taken from test_shallow_water
	1710	"""
	1711	quantity = Quantity(self.mesh4)
	1712	quantity.domain.beta_w = 0.9
	1713
	1714	#Test centroids
	1715	quantity.set_values([2.,4.,8.,2.], location = 'centroids')
	1716	assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
	1717
	1718
	1719	#Extrapolate
	1720	quantity.extrapolate_second_order()
	1721
	1722	assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
	1723
	1724	#Limit
	1725	quantity.limit()
	1726
	1727	# limited value for beta_w = 0.9
	1728
	1729	assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
	1730	# limited values for beta_w = 0.5
	1731	#assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
	1732
	1733
	1734	#Assert that quantities are conserved
	1735	from Numeric import sum
	1736	for k in range(quantity.centroid_values.shape[0]):
	1737	assert allclose (quantity.centroid_values[k],
	1738	sum(quantity.vertex_values[k,:])/3)
	1739
	1740
	1741
	1742
	1743
	1744	def test_distribute_first_order(self):
	1745	quantity = Quantity(self.mesh4)
	1746
	1747	#Test centroids
	1748	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
	1749	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
	1750
	1751
	1752	#Extrapolate from centroid to vertices and edges
	1753	quantity.extrapolate_first_order()
	1754
	1755	#Interpolate
	1756	#quantity.interpolate_from_vertices_to_edges()
	1757
	1758	assert allclose(quantity.vertex_values,
	1759	[[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
	1760	assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
	1761	[3,3,3], [4, 4, 4]])
	1762
	1763
	1764	def test_interpolate_from_vertices_to_edges(self):
	1765	quantity = Quantity(self.mesh4)
	1766
	1767	quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float)
	1768
	1769	quantity.interpolate_from_vertices_to_edges()
	1770
	1771	assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
	1772	[3., 2.5, 1.5],
	1773	[3.5, 4.5, 3.],
	1774	[2.5, 3.5, 2]])
	1775
	1776
	1777	def test_interpolate_from_edges_to_vertices(self):
	1778	quantity = Quantity(self.mesh4)
	1779
	1780	quantity.edge_values = array([[1., 1.5, 0.5],
	1781	[3., 2.5, 1.5],
	1782	[3.5, 4.5, 3.],
	1783	[2.5, 3.5, 2]],Float)
	1784
	1785	quantity.interpolate_from_edges_to_vertices()
	1786
	1787	assert allclose(quantity.vertex_values,
	1788	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
	1789
	1790
	1791
	1792	def test_distribute_second_order(self):
	1793	quantity = Quantity(self.mesh4)
	1794
	1795	#Test centroids
	1796	quantity.set_values([2.,4.,8.,2.], location = 'centroids')
	1797	assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
	1798
	1799
	1800	#Extrapolate
	1801	quantity.extrapolate_second_order()
	1802
	1803	assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
	1804
	1805
	1806	def test_update_explicit(self):
	1807	quantity = Quantity(self.mesh4)
	1808
	1809	#Test centroids
	1810	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
	1811	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
	1812
	1813	#Set explicit_update
	1814	quantity.explicit_update = array( [1.,1.,1.,1.] )
	1815
	1816	#Update with given timestep
	1817	quantity.update(0.1)
	1818
	1819	x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
	1820	assert allclose( quantity.centroid_values, x)
	1821
	1822	def test_update_semi_implicit(self):
	1823	quantity = Quantity(self.mesh4)
	1824
	1825	#Test centroids
	1826	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
	1827	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
	1828
	1829	#Set semi implicit update
	1830	quantity.semi_implicit_update = array([1.,1.,1.,1.])
	1831
	1832	#Update with given timestep
	1833	timestep = 0.1
	1834	quantity.update(timestep)
	1835
	1836	sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
	1837	denom = ones(4, Float)-timestep*sem
	1838
	1839	x = array([1, 2, 3, 4])/denom
	1840	assert allclose( quantity.centroid_values, x)
	1841
	1842
	1843	def test_both_updates(self):
	1844	quantity = Quantity(self.mesh4)
	1845
	1846	#Test centroids
	1847	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
	1848	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
	1849
	1850	#Set explicit_update
	1851	quantity.explicit_update = array( [4.,3.,2.,1.] )
	1852
	1853	#Set semi implicit update
	1854	quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
	1855
	1856	#Update with given timestep
	1857	timestep = 0.1
	1858	quantity.update(0.1)
	1859
	1860	sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
	1861	denom = ones(4, Float)-timestep*sem
	1862
	1863	x = array([1., 2., 3., 4.])
	1864	x /= denom
	1865	x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
	1866
	1867	assert allclose( quantity.centroid_values, x)
	1868
	1869
	1870
	1871
	1872	#Test smoothing
	1873	def test_smoothing(self):
	1874
	1875	from mesh_factory import rectangular
	1876	from shallow_water import Domain, Transmissive_boundary
	1877	from Numeric import zeros, Float
	1878	from anuga.utilities.numerical_tools import mean
	1879
	1880	#Create basic mesh
	1881	points, vertices, boundary = rectangular(2, 2)
	1882
	1883	#Create shallow water domain
	1884	domain = Domain(points, vertices, boundary)
	1885	domain.default_order=2
	1886	domain.reduction = mean
	1887
	1888
	1889	#Set some field values
	1890	domain.set_quantity('elevation', lambda x,y: x)
	1891	domain.set_quantity('friction', 0.03)
	1892
	1893
	1894	######################
	1895	# Boundary conditions
	1896	B = Transmissive_boundary(domain)
	1897	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
	1898
	1899
	1900	######################
	1901	#Initial condition - with jumps
	1902
	1903	bed = domain.quantities['elevation'].vertex_values
	1904	stage = zeros(bed.shape, Float)
	1905
	1906	h = 0.03
	1907	for i in range(stage.shape[0]):
	1908	if i % 2 == 0:
	1909	stage[i,:] = bed[i,:] + h
	1910	else:
	1911	stage[i,:] = bed[i,:]
	1912
	1913	domain.set_quantity('stage', stage)
	1914
	1915	stage = domain.quantities['stage']
	1916
	1917	#Get smoothed stage
	1918	A, V = stage.get_vertex_values(xy=False, smooth=True)
	1919	Q = stage.vertex_values
	1920
	1921
	1922	assert A.shape[0] == 9
	1923	assert V.shape[0] == 8
	1924	assert V.shape[1] == 3
	1925
	1926	#First four points
	1927	assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
	1928	assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
	1929	assert allclose(A[2], Q[3,0])
	1930	assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
	1931
	1932	#Center point
	1933	assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
	1934	Q[5,0] + Q[6,2] + Q[7,1])/6)
	1935
	1936
	1937	#Check V
	1938	assert allclose(V[0,:], [3,4,0])
	1939	assert allclose(V[1,:], [1,0,4])
	1940	assert allclose(V[2,:], [4,5,1])
	1941	assert allclose(V[3,:], [2,1,5])
	1942	assert allclose(V[4,:], [6,7,3])
	1943	assert allclose(V[5,:], [4,3,7])
	1944	assert allclose(V[6,:], [7,8,4])
	1945	assert allclose(V[7,:], [5,4,8])
	1946
	1947	#Get smoothed stage with XY
	1948	X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
	1949
	1950	assert allclose(A, A1)
	1951	assert allclose(V, V1)
	1952
	1953	#Check XY
	1954	assert allclose(X[4], 0.5)
	1955	assert allclose(Y[4], 0.5)
	1956
	1957	assert allclose(X[7], 1.0)
	1958	assert allclose(Y[7], 0.5)
	1959
	1960
	1961
	1962
	1963	def test_vertex_values_no_smoothing(self):
	1964
	1965	from mesh_factory import rectangular
	1966	from shallow_water import Domain, Transmissive_boundary
	1967	from Numeric import zeros, Float
	1968	from anuga.utilities.numerical_tools import mean
	1969
	1970
	1971	#Create basic mesh
	1972	points, vertices, boundary = rectangular(2, 2)
	1973
	1974	#Create shallow water domain
	1975	domain = Domain(points, vertices, boundary)
	1976	domain.default_order=2
	1977	domain.reduction = mean
	1978
	1979
	1980	#Set some field values
	1981	domain.set_quantity('elevation', lambda x,y: x)
	1982	domain.set_quantity('friction', 0.03)
	1983
	1984
	1985	######################
	1986	#Initial condition - with jumps
	1987
	1988	bed = domain.quantities['elevation'].vertex_values
	1989	stage = zeros(bed.shape, Float)
	1990
	1991	h = 0.03
	1992	for i in range(stage.shape[0]):
	1993	if i % 2 == 0:
	1994	stage[i,:] = bed[i,:] + h
	1995	else:
	1996	stage[i,:] = bed[i,:]
	1997
	1998	domain.set_quantity('stage', stage)
	1999
	2000	#Get stage
	2001	stage = domain.quantities['stage']
	2002	A, V = stage.get_vertex_values(xy=False, smooth=False)
	2003	Q = stage.vertex_values.flat
	2004
	2005	for k in range(8):
	2006	assert allclose(A[k], Q[k])
	2007
	2008
	2009	for k in range(8):
	2010	assert V[k, 0] == 3*k
	2011	assert V[k, 1] == 3*k+1
	2012	assert V[k, 2] == 3*k+2
	2013
	2014
	2015
	2016	X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
	2017
	2018
	2019	assert allclose(A, A1)
	2020	assert allclose(V, V1)
	2021
	2022	#Check XY
	2023	assert allclose(X[1], 0.5)
	2024	assert allclose(Y[1], 0.5)
	2025	assert allclose(X[4], 0.0)
	2026	assert allclose(Y[4], 0.0)
	2027	assert allclose(X[12], 1.0)
	2028	assert allclose(Y[12], 0.0)
	2029
	2030
	2031
	2032	def set_array_values_by_index(self):
	2033
	2034	from mesh_factory import rectangular
	2035	from shallow_water import Domain
	2036	from Numeric import zeros, Float
	2037
	2038	#Create basic mesh
	2039	points, vertices, boundary = rectangular(1, 1)
	2040
	2041	#Create shallow water domain
	2042	domain = Domain(points, vertices, boundary)
	2043	#print "domain.number_of_elements ",domain.number_of_elements
	2044	quantity = Quantity(domain,[[1,1,1],[2,2,2]])
	2045	value = [7]
	2046	indices = [1]
	2047	quantity.set_array_values_by_index(value,
	2048	location = 'centroids',
	2049	indices = indices)
	2050	#print "quantity.centroid_values",quantity.centroid_values
	2051
	2052	assert allclose(quantity.centroid_values, [1,7])
	2053
	2054	quantity.set_array_values([15,20,25], indices = indices)
	2055	assert allclose(quantity.centroid_values, [1,20])
	2056
	2057	quantity.set_array_values([15,20,25], indices = indices)
	2058	assert allclose(quantity.centroid_values, [1,20])
	2059
	2060	def test_setting_some_vertex_values(self):
	2061	"""
	2062	set values based on triangle lists.
	2063	"""
	2064	from mesh_factory import rectangular
	2065	from shallow_water import Domain
	2066	from Numeric import zeros, Float
	2067
	2068	#Create basic mesh
	2069	points, vertices, boundary = rectangular(1, 3)
	2070	#print "vertices",vertices
	2071	#Create shallow water domain
	2072	domain = Domain(points, vertices, boundary)
	2073	#print "domain.number_of_elements ",domain.number_of_elements
	2074	quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
	2075	[4,4,4],[5,5,5],[6,6,6]])
	2076
	2077
	2078	# Check that constants work
	2079	value = 7
	2080	indices = [1]
	2081	quantity.set_values(value,
	2082	location = 'centroids',
	2083	indices = indices)
	2084	#print "quantity.centroid_values",quantity.centroid_values
	2085	assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
	2086
	2087	value = [7]
	2088	indices = [1]
	2089	quantity.set_values(value,
	2090	location = 'centroids',
	2091	indices = indices)
	2092	#print "quantity.centroid_values",quantity.centroid_values
	2093	assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
	2094
	2095	value = [[15,20,25]]
	2096	quantity.set_values(value, indices = indices)
	2097	#print "1 quantity.vertex_values",quantity.vertex_values
	2098	assert allclose(quantity.vertex_values[1], value[0])
	2099
	2100
	2101	#print "quantity",quantity.vertex_values
	2102	values = [10,100,50]
	2103	quantity.set_values(values, indices = [0,1,5], location = 'centroids')
	2104	#print "2 quantity.vertex_values",quantity.vertex_values
	2105	assert allclose(quantity.vertex_values[0], [10,10,10])
	2106	assert allclose(quantity.vertex_values[5], [50,50,50])
	2107	#quantity.interpolate()
	2108	#print "quantity.centroid_values",quantity.centroid_values
	2109	assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
	2110
	2111
	2112	quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
	2113	[4,4,4],[5,5,5],[6,6,6]])
	2114	values = [10,100,50]
	2115	#this will be per unique vertex, indexing the vertices
	2116	#print "quantity.vertex_values",quantity.vertex_values
	2117	quantity.set_values(values, indices = [0,1,5])
	2118	#print "quantity.vertex_values",quantity.vertex_values
	2119	assert allclose(quantity.vertex_values[0], [1,50,10])
	2120	assert allclose(quantity.vertex_values[5], [6,6,6])
	2121	assert allclose(quantity.vertex_values[1], [100,10,50])
	2122
	2123	quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
	2124	[4,4,4],[5,5,5],[6,6,6]])
	2125	values = [[31,30,29],[400,400,400],[1000,999,998]]
	2126	quantity.set_values(values, indices = [3,3,5])
	2127	quantity.interpolate()
	2128	assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
	2129
	2130	values = [[1,1,1],[2,2,2],[3,3,3],
	2131	[4,4,4],[5,5,5],[6,6,6]]
	2132	quantity.set_values(values)
	2133
	2134	# testing the standard set values by vertex
	2135	# indexed by vertex_id in general_mesh.coordinates
	2136	values = [0,1,2,3,4,5,6,7]
	2137
	2138	quantity.set_values(values)
	2139	#print "1 quantity.vertex_values",quantity.vertex_values
	2140	assert allclose(quantity.vertex_values,[[ 4., 5., 0.],
	2141	[ 1., 0., 5.],
	2142	[ 5., 6., 1.],
	2143	[ 2., 1., 6.],
	2144	[ 6., 7., 2.],
	2145	[ 3., 2., 7.]])
	2146
	2147	def test_setting_unique_vertex_values(self):
	2148	"""
	2149	set values based on unique_vertex lists.
	2150	"""
	2151	from mesh_factory import rectangular
	2152	from shallow_water import Domain
	2153	from Numeric import zeros, Float
	2154
	2155	#Create basic mesh
	2156	points, vertices, boundary = rectangular(1, 3)
	2157	#print "vertices",vertices
	2158	#Create shallow water domain
	2159	domain = Domain(points, vertices, boundary)
	2160	#print "domain.number_of_elements ",domain.number_of_elements
	2161	quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
	2162	[4,4,4],[5,5,5]])
	2163	value = 7
	2164	indices = [1,5]
	2165	quantity.set_values(value,
	2166	location = 'unique vertices',
	2167	indices = indices)
	2168	#print "quantity.centroid_values",quantity.centroid_values
	2169	assert allclose(quantity.vertex_values[0], [0,7,0])
	2170	assert allclose(quantity.vertex_values[1], [7,1,7])
	2171	assert allclose(quantity.vertex_values[2], [7,2,7])
	2172
	2173
	2174	def test_get_values(self):
	2175	"""
	2176	get values based on triangle lists.
	2177	"""
	2178	from mesh_factory import rectangular
	2179	from shallow_water import Domain
	2180	from Numeric import zeros, Float
	2181
	2182	#Create basic mesh
	2183	points, vertices, boundary = rectangular(1, 3)
	2184
	2185	#print "points",points
	2186	#print "vertices",vertices
	2187	#print "boundary",boundary
	2188
	2189	#Create shallow water domain
	2190	domain = Domain(points, vertices, boundary)
	2191	#print "domain.number_of_elements ",domain.number_of_elements
	2192	quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
	2193	[4,4,4],[5,5,5]])
	2194
	2195	#print "quantity.get_values(location = 'unique vertices')", \
	2196	# quantity.get_values(location = 'unique vertices')
	2197
	2198	#print "quantity.get_values(location = 'unique vertices')", \
	2199	# quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
	2200	# location = 'unique vertices')
	2201
	2202	answer = [0.5,2,4,5,0,1,3,4.5]
	2203	assert allclose(answer,
	2204	quantity.get_values(location = 'unique vertices'))
	2205
	2206	indices = [0,5,3]
	2207	answer = [0.5,1,5]
	2208	assert allclose(answer,
	2209	quantity.get_values(indices=indices, \
	2210	location = 'unique vertices'))
	2211	#print "quantity.centroid_values",quantity.centroid_values
	2212	#print "quantity.get_values(location = 'centroids') ",\
	2213	# quantity.get_values(location = 'centroids')
	2214
	2215
	2216
	2217
	2218	def test_get_values_2(self):
	2219	"""Different mesh (working with domain object) - also check centroids.
	2220	"""
	2221
	2222
	2223	a = [0.0, 0.0]
	2224	b = [0.0, 2.0]
	2225	c = [2.0,0.0]
	2226	d = [0.0, 4.0]
	2227	e = [2.0, 2.0]
	2228	f = [4.0,0.0]
	2229
	2230	points = [a, b, c, d, e, f]
	2231	#bac, bce, ecf, dbe
	2232	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	2233
	2234	domain = Domain(points, vertices)
	2235
	2236	quantity = Quantity(domain)
	2237	quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
	2238
	2239	assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
	2240	assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
	2241
	2242
	2243	assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
	2244	[4,2,6],
	2245	[6,2,4],
	2246	[8,4,6]])
	2247
	2248	assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
	2249	[8,4,6]])
	2250
	2251
	2252	assert allclose(quantity.get_values(location='edges'), [[1,3,2],
	2253	[4,5,3],
	2254	[3,5,4],
	2255	[5,7,6]])
	2256	assert allclose(quantity.get_values(location='edges', indices=[1,3]),
	2257	[[4,5,3],
	2258	[5,7,6]])
	2259
	2260	# Check averaging over vertices
	2261	#a: 0
	2262	#b: (4+4+4)/3
	2263	#c: (2+2+2)/3
	2264	#d: 8
	2265	#e: (6+6+6)/3
	2266	#f: 4
	2267	assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4])
	2268
	2269
	2270
	2271
	2272
	2273
	2274	def test_get_interpolated_values(self):
	2275
	2276	from mesh_factory import rectangular
	2277	from shallow_water import Domain
	2278	from Numeric import zeros, Float
	2279
	2280	#Create basic mesh
	2281	points, vertices, boundary = rectangular(1, 3)
	2282	domain = Domain(points, vertices, boundary)
	2283
	2284	#Constant values
	2285	quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
	2286	[4,4,4],[5,5,5]])
	2287
	2288
	2289
	2290	# Get interpolated values at centroids
	2291	interpolation_points = domain.get_centroid_coordinates()
	2292	answer = quantity.get_values(location='centroids')
	2293
	2294
	2295	#print quantity.get_values(points=interpolation_points)
	2296	assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
	2297
	2298
	2299	#Arbitrary values
	2300	quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
	2301	[1,4,-9],[2,5,0]])
	2302
	2303
	2304	# Get interpolated values at centroids
	2305	interpolation_points = domain.get_centroid_coordinates()
	2306	answer = quantity.get_values(location='centroids')
	2307	#print answer
	2308	#print quantity.get_values(points=interpolation_points)
	2309	assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
	2310
	2311
	2312	#FIXME TODO
	2313	#indices = [0,5,3]
	2314	#answer = [0.5,1,5]
	2315	#assert allclose(answer,
	2316	# quantity.get_values(indices=indices, \
	2317	# location = 'unique vertices'))
	2318
	2319
	2320
	2321
	2322	def test_get_interpolated_values_2(self):
	2323	a = [0.0, 0.0]
	2324	b = [0.0, 2.0]
	2325	c = [2.0,0.0]
	2326	d = [0.0, 4.0]
	2327	e = [2.0, 2.0]
	2328	f = [4.0,0.0]
	2329
	2330	points = [a, b, c, d, e, f]
	2331	#bac, bce, ecf, dbe
	2332	vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	2333
	2334	domain = Domain(points, vertices)
	2335
	2336	quantity = Quantity(domain)
	2337	quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
	2338
	2339	#First pick one point
	2340	x, y = 2.0/3, 8.0/3
	2341	v = quantity.get_values(interpolation_points = [[x,y]])
	2342	assert allclose(v, 6)
	2343
	2344	# Then another to test that algorithm won't blindly
	2345	# reuse interpolation matrix
	2346	x, y = 4.0/3, 4.0/3
	2347	v = quantity.get_values(interpolation_points = [[x,y]])
	2348	assert allclose(v, 4)
	2349
	2350
	2351
	2352
	2353	def test_getting_some_vertex_values(self):
	2354	"""
	2355	get values based on triangle lists.
	2356	"""
	2357	from mesh_factory import rectangular
	2358	from shallow_water import Domain
	2359	from Numeric import zeros, Float
	2360
	2361	#Create basic mesh
	2362	points, vertices, boundary = rectangular(1, 3)
	2363
	2364	#print "points",points
	2365	#print "vertices",vertices
	2366	#print "boundary",boundary
	2367
	2368	#Create shallow water domain
	2369	domain = Domain(points, vertices, boundary)
	2370	#print "domain.number_of_elements ",domain.number_of_elements
	2371	quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
	2372	[4,4,4],[5,5,5],[6,6,6]])
	2373	value = [7]
	2374	indices = [1]
	2375	quantity.set_values(value,
	2376	location = 'centroids',
	2377	indices = indices)
	2378	#print "quantity.centroid_values",quantity.centroid_values
	2379	#print "quantity.get_values(location = 'centroids') ",\
	2380	# quantity.get_values(location = 'centroids')
	2381	assert allclose(quantity.centroid_values,
	2382	quantity.get_values(location = 'centroids'))
	2383
	2384
	2385	value = [[15,20,25]]
	2386	quantity.set_values(value, indices = indices)
	2387	#print "1 quantity.vertex_values",quantity.vertex_values
	2388	assert allclose(quantity.vertex_values, quantity.get_values())
	2389
	2390	assert allclose(quantity.edge_values,
	2391	quantity.get_values(location = 'edges'))
	2392
	2393	# get a subset of elements
	2394	subset = quantity.get_values(location='centroids', indices=[0,5])
	2395	answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
	2396	assert allclose(subset, answer)
	2397
	2398
	2399	subset = quantity.get_values(location='edges', indices=[0,5])
	2400	answer = [quantity.edge_values[0],quantity.edge_values[5]]
	2401	#print "subset",subset
	2402	#print "answer",answer
	2403	assert allclose(subset, answer)
	2404
	2405	subset = quantity.get_values( indices=[1,5])
	2406	answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
	2407	#print "subset",subset
	2408	#print "answer",answer
	2409	assert allclose(subset, answer)
	2410
	2411	def test_smooth_vertex_values(self):
	2412	"""
	2413	get values based on triangle lists.
	2414	"""
	2415	from mesh_factory import rectangular
	2416	from shallow_water import Domain
	2417	from Numeric import zeros, Float
	2418
	2419	#Create basic mesh
	2420	points, vertices, boundary = rectangular(2, 2)
	2421
	2422	#print "points",points
	2423	#print "vertices",vertices
	2424	#print "boundary",boundary
	2425
	2426	#Create shallow water domain
	2427	domain = Domain(points, vertices, boundary)
	2428	#print "domain.number_of_elements ",domain.number_of_elements
	2429	quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
	2430	[4,4,4],[5,5,5],[6,6,6],[7,7,7]])
	2431
	2432	#print "quantity.get_values(location = 'unique vertices')", \
	2433	# quantity.get_values(location = 'unique vertices')
	2434
	2435	#print "quantity.get_values(location = 'unique vertices')", \
	2436	# quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
	2437	# location = 'unique vertices')
	2438
	2439	#print quantity.get_values(location = 'unique vertices')
	2440	#print quantity.domain.number_of_triangles_per_node
	2441	#print quantity.vertex_values
	2442
	2443	#answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5]
	2444	#assert allclose(answer,
	2445	# quantity.get_values(location = 'unique vertices'))
	2446
	2447	quantity.smooth_vertex_values()
	2448
	2449	#print quantity.vertex_values
	2450
	2451
	2452	answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4],
	2453	[4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]]
	2454
	2455	assert allclose(answer_vertex_values,
	2456	quantity.vertex_values)
	2457	#print "quantity.centroid_values",quantity.centroid_values
	2458	#print "quantity.get_values(location = 'centroids') ",\
	2459	# quantity.get_values(location = 'centroids')
	2460
	2461
	2462
	2463	#-------------------------------------------------------------
	2464	if __name__ == "__main__":
	2465	suite = unittest.makeSuite(Test_Quantity, 'test')
	2466	#suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
	2467
	2468	#suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
	2469	#print "restricted test"
	2470	#suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
	2471	runner = unittest.TextTestRunner()
	2472	runner.run(suite)

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

Download in other formats: