Context Navigation

source: anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py @ 5517

Last change on this file since 5517 was 5517, checked in by ole, 15 years ago
Work on test_quantity towards general ability to use a polygon.
File size: 81.6 KB

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

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

Download in other formats: