Context Navigation

source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py @ 8063

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

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

Download in other formats: