Context Navigation

test_quantity.py @ 6883

Last change on this file since 6883 was 6410, checked in by rwilson, 16 years ago
numpy changes.
File size: 85.8 KB

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

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

Context Navigation

source: branches/numpy/anuga/abstract_2d_finite_volumes/test_quantity.py @ 6883

Download in other formats: