Context Navigation

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

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

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

Download in other formats: