Context Navigation

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

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

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

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

Download in other formats: