Context Navigation

source: inundation/pyvolution/test_quantity.py @ 1869

Last change on this file since 1869 was 1754, checked in by ole, 19 years ago
Implemented operator overloading for quantities Implemented set_values with quantity as input + testing
File size: 41.4 KB

Line
1	#!/usr/bin/env python
2
3	import unittest
4	from math import sqrt, pi
5
6
7	from quantity import *
8	from config import epsilon
9	from Numeric import allclose, array, ones, Float
10
11
12	#Aux for least_squares example
13	def linear_function(point):
14	point = array(point)
15	return point[:,0]+point[:,1]
16
17
18	class Test_Quantity(unittest.TestCase):
19	def setUp(self):
20	from domain import Domain
21
22	a = [0.0, 0.0]
23	b = [0.0, 2.0]
24	c = [2.0, 0.0]
25	d = [0.0, 4.0]
26	e = [2.0, 2.0]
27	f = [4.0, 0.0]
28
29	points = [a, b, c, d, e, f]
30
31	#bac, bce, ecf, dbe
32	elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
33
34	self.mesh1 = Domain(points[:3], [elements[0]])
35	self.mesh1.check_integrity()
36
37	self.mesh4 = Domain(points, elements)
38	self.mesh4.check_integrity()
39
40	def tearDown(self):
41	pass
42	#print " Tearing down"
43
44
45	def test_creation(self):
46
47	quantity = Quantity(self.mesh1, [[1,2,3]])
48	assert allclose(quantity.vertex_values, [[1.,2.,3.]])
49
50	try:
51	quantity = Quantity()
52	except:
53	pass
54	else:
55	raise 'Should have raised empty quantity exception'
56
57
58	try:
59	quantity = Quantity([1,2,3])
60	except AssertionError:
61	pass
62	except:
63	raise 'Should have raised "mising mesh object" error'
64
65
66	def test_creation_zeros(self):
67
68	quantity = Quantity(self.mesh1)
69	assert allclose(quantity.vertex_values, [[0.,0.,0.]])
70
71
72	quantity = Quantity(self.mesh4)
73	assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
74	[0.,0.,0.], [0.,0.,0.]])
75
76
77	def test_interpolation(self):
78	quantity = Quantity(self.mesh1, [[1,2,3]])
79	assert allclose(quantity.centroid_values, [2.0]) #Centroid
80
81	assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
82
83
84	def test_interpolation2(self):
85	quantity = Conserved_quantity(self.mesh4,
86	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
87	assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
88
89
90	quantity.extrapolate_second_order()
91
92	#print quantity.vertex_values
93	#assert allclose(quantity.vertex_values, [[2., 2., 2.],
94	# [3.+2./3, 6.+2./3, 4.+2./3],
95	# [7.5, 0.5, 1.],
96	# [-5, -2.5, 7.5]])
97
98	assert allclose(quantity.vertex_values[1,:],[3.+2./3, 6.+2./3, 4.+2./3])
99	#FIXME: Work out the others
100
101
102	#print quantity.edge_values
103	assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
104	[5., 5., 5.],
105	[4.5, 4.5, 0.],
106	[3.0, -1.5, -1.5]])
107
108	def test_boundary_allocation(self):
109	quantity = Conserved_quantity(self.mesh4,
110	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
111
112	assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
113
114
115	def test_set_values(self):
116	quantity = Quantity(self.mesh4)
117
118
119	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
120	location = 'vertices')
121	assert allclose(quantity.vertex_values,
122	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
123	assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
124	assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
125	[5., 5., 5.],
126	[4.5, 4.5, 0.],
127	[3.0, -1.5, -1.5]])
128
129
130	#Test default
131	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
132	assert allclose(quantity.vertex_values,
133	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
134	assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
135	assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
136	[5., 5., 5.],
137	[4.5, 4.5, 0.],
138	[3.0, -1.5, -1.5]])
139
140	#Test centroids
141	quantity.set_values([1,2,3,4], location = 'centroids')
142	assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
143
144	#Test edges
145	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
146	location = 'edges')
147	assert allclose(quantity.edge_values,
148	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
149
150	#Test exceptions
151	try:
152	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
153	location = 'bas kamel tuba')
154	except:
155	pass
156
157
158	try:
159	quantity.set_values([[1,2,3], [0,0,9]])
160	except AssertionError:
161	pass
162	except:
163	raise 'should have raised Assertionerror'
164
165
166
167	def test_set_values_const(self):
168	quantity = Quantity(self.mesh4)
169
170	quantity.set_values(1.0, location = 'vertices')
171	assert allclose(quantity.vertex_values,
172	[[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
173
174	assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
175	assert allclose(quantity.edge_values, [[1, 1, 1],
176	[1, 1, 1],
177	[1, 1, 1],
178	[1, 1, 1]])
179
180
181	quantity.set_values(2.0, location = 'centroids')
182	assert allclose(quantity.centroid_values, [2, 2, 2, 2])
183
184	quantity.set_values(3.0, location = 'edges')
185	assert allclose(quantity.edge_values, [[3, 3, 3],
186	[3, 3, 3],
187	[3, 3, 3],
188	[3, 3, 3]])
189
190
191	def test_set_values_func(self):
192	quantity = Quantity(self.mesh4)
193
194	def f(x, y):
195	return x+y
196
197	quantity.set_values(f, location = 'vertices')
198	#print "quantity.vertex_values",quantity.vertex_values
199	assert allclose(quantity.vertex_values,
200	[[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
201	assert allclose(quantity.centroid_values,
202	[4.0/3, 8.0/3, 10.0/3, 10.0/3])
203	assert allclose(quantity.edge_values,
204	[[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
205
206
207	quantity.set_values(f, location = 'centroids')
208	assert allclose(quantity.centroid_values,
209	[4.0/3, 8.0/3, 10.0/3, 10.0/3])
210
211
212	def test_integral(self):
213	quantity = Quantity(self.mesh4)
214
215	#Try constants first
216	const = 5
217	quantity.set_values(const, location = 'vertices')
218	#print 'Q', quantity.get_integral()
219
220	assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)
221
222	#Try with a linear function
223	def f(x, y):
224	return x+y
225
226	quantity.set_values(f, location = 'vertices')
227
228
229	ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
230
231	assert allclose (quantity.get_integral(), ref_integral)
232
233
234
235	def test_set_vertex_values(self):
236	quantity = Quantity(self.mesh4)
237	quantity.set_vertex_values([0,1,2,3,4,5])
238
239	assert allclose(quantity.vertex_values,
240	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
241	assert allclose(quantity.centroid_values,
242	[1., 7./3, 11./3, 8./3]) #Centroid
243	assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
244	[3., 2.5, 1.5],
245	[3.5, 4.5, 3.],
246	[2.5, 3.5, 2]])
247
248
249	def test_set_vertex_values_subset(self):
250	quantity = Quantity(self.mesh4)
251	quantity.set_vertex_values([0,1,2,3,4,5])
252	quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
253
254	assert allclose(quantity.vertex_values,
255	[[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
256
257
258	def test_set_vertex_values_using_general_interface(self):
259	quantity = Quantity(self.mesh4)
260
261
262	quantity.set_values([0,1,2,3,4,5])
263
264
265	assert allclose(quantity.vertex_values,
266	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
267
268	#Centroid
269	assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
270
271	assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
272	[3., 2.5, 1.5],
273	[3.5, 4.5, 3.],
274	[2.5, 3.5, 2]])
275
276
277
278
279
280	def test_set_values_using_least_squares(self):
281
282
283	quantity = Quantity(self.mesh4)
284
285	#Get (enough) datapoints
286	data_points = [[ 0.66666667, 0.66666667],
287	[ 1.33333333, 1.33333333],
288	[ 2.66666667, 0.66666667],
289	[ 0.66666667, 2.66666667],
290	[ 0.0, 1.0],
291	[ 0.0, 3.0],
292	[ 1.0, 0.0],
293	[ 1.0, 1.0],
294	[ 1.0, 2.0],
295	[ 1.0, 3.0],
296	[ 2.0, 1.0],
297	[ 3.0, 0.0],
298	[ 3.0, 1.0]]
299
300	z = linear_function(data_points)
301
302	#Obsoleted bit
303	#interp = Interpolation(quantity.domain.coordinates,
304	# quantity.domain.triangles,
305	# data_points, alpha=0.0)
306	#
307	#answer = linear_function(quantity.domain.coordinates)
308	#
309	#f = interp.fit(z)
310	#
311	#print "f",f
312	#print "answer",answer
313	#assert allclose(f, answer)
314
315
316	#Use built-in least squares fit
317	quantity.set_values(points = data_points, values = z, alpha = 0)
318	answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True))
319	#print quantity.vertex_values, answer
320	assert allclose(quantity.vertex_values.flat, answer)
321
322
323	#Now try by setting the same values directly
324	from least_squares import fit_to_mesh
325	vertex_attributes = fit_to_mesh(quantity.domain.coordinates,
326	quantity.domain.triangles,
327	data_points,
328	z,
329	alpha = 0,
330	verbose=False)
331
332	#print vertex_attributes
333	quantity.set_values(vertex_attributes)
334	assert allclose(quantity.vertex_values.flat, answer)
335
336
337	def test_set_values_from_file(self):
338	quantity = Quantity(self.mesh4)
339
340	#Get (enough) datapoints
341	data_points = [[ 0.66666667, 0.66666667],
342	[ 1.33333333, 1.33333333],
343	[ 2.66666667, 0.66666667],
344	[ 0.66666667, 2.66666667],
345	[ 0.0, 1.0],
346	[ 0.0, 3.0],
347	[ 1.0, 0.0],
348	[ 1.0, 1.0],
349	[ 1.0, 2.0],
350	[ 1.0, 3.0],
351	[ 2.0, 1.0],
352	[ 3.0, 0.0],
353	[ 3.0, 1.0]]
354
355	z = linear_function(data_points)
356
357
358	#Create pts file
359	from data_manager import write_ptsfile
360	ptsfile = 'testptsfile.pts'
361	att = 'spam_and_eggs'
362	write_ptsfile(ptsfile, data_points, z, attribute_name = att)
363
364	#Check that values can be set from file
365	quantity.set_values(filename = ptsfile,
366	attribute_name = att, alpha = 0)
367	answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True))
368	#print quantity.vertex_values, answer
369	assert allclose(quantity.vertex_values.flat, answer)
370
371
372	#Check that values can be set from file using default attribute
373	quantity.set_values(filename = ptsfile, alpha = 0)
374	assert allclose(quantity.vertex_values.flat, answer)
375
376	#Cleanup
377	import os
378	os.remove(ptsfile)
379
380
381
382	def test_set_values_from_quantity(self):
383
384	quantity1 = Quantity(self.mesh4)
385	quantity1.set_vertex_values([0,1,2,3,4,5])
386
387	assert allclose(quantity1.vertex_values,
388	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
389
390
391	quantity2 = Quantity(self.mesh4)
392	quantity2.set_values(quantity = quantity1)
393	assert allclose(quantity2.vertex_values,
394	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
395
396	quantity2.set_values(quantity = 2*quantity1)
397	assert allclose(quantity2.vertex_values,
398	[[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
399
400	quantity2.set_values(quantity = 2*quantity1 + 3)
401	assert allclose(quantity2.vertex_values,
402	[[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
403
404
405	#Check detection of quantity as first orgument
406	quantity2.set_values(2*quantity1 + 3)
407	assert allclose(quantity2.vertex_values,
408	[[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
409
410
411
412
413
414	def test_overloading(self):
415
416	quantity1 = Quantity(self.mesh4)
417	quantity1.set_vertex_values([0,1,2,3,4,5])
418
419	assert allclose(quantity1.vertex_values,
420	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
421
422
423	quantity2 = Quantity(self.mesh4)
424	quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
425	location = 'vertices')
426
427
428
429	quantity3 = Quantity(self.mesh4)
430	quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]],
431	location = 'vertices')
432
433
434	#Negation
435	Q = -quantity1
436	assert allclose(Q.vertex_values, -quantity1.vertex_values)
437	assert allclose(Q.centroid_values, -quantity1.centroid_values)
438	assert allclose(Q.edge_values, -quantity1.edge_values)
439
440	#Addition
441	Q = quantity1 + 7
442	assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
443	assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
444	assert allclose(Q.edge_values, quantity1.edge_values + 7)
445
446	Q = 7 + quantity1
447	assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
448	assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
449	assert allclose(Q.edge_values, quantity1.edge_values + 7)
450
451	Q = quantity1 + quantity2
452	assert allclose(Q.vertex_values,
453	quantity1.vertex_values + quantity2.vertex_values)
454	assert allclose(Q.centroid_values,
455	quantity1.centroid_values + quantity2.centroid_values)
456	assert allclose(Q.edge_values,
457	quantity1.edge_values + quantity2.edge_values)
458
459
460	Q = quantity1 + quantity2 - 3
461	assert allclose(Q.vertex_values,
462	quantity1.vertex_values + quantity2.vertex_values - 3)
463
464	Q = quantity1 - quantity2
465	assert allclose(Q.vertex_values,
466	quantity1.vertex_values - quantity2.vertex_values)
467
468	#Scaling
469	Q = quantity1*3
470	assert allclose(Q.vertex_values, quantity1.vertex_values*3)
471	assert allclose(Q.centroid_values, quantity1.centroid_values*3)
472	assert allclose(Q.edge_values, quantity1.edge_values*3)
473	Q = 3*quantity1
474	assert allclose(Q.vertex_values, quantity1.vertex_values*3)
475
476	#Multiplication
477	Q = quantity1 * quantity2
478	#print Q.vertex_values
479	#print Q.centroid_values
480	#print quantity1.centroid_values
481	#print quantity2.centroid_values
482
483	assert allclose(Q.vertex_values,
484	quantity1.vertex_values * quantity2.vertex_values)
485
486	#Linear combinations
487	Q = 4*quantity1 + 2
488	assert allclose(Q.vertex_values,
489	4*quantity1.vertex_values + 2)
490
491	Q = quantity1*quantity2 + 2
492	assert allclose(Q.vertex_values,
493	quantity1.vertex_values * quantity2.vertex_values + 2)
494
495	Q = quantity1*quantity2 + quantity3
496	assert allclose(Q.vertex_values,
497	quantity1.vertex_values * quantity2.vertex_values +
498	quantity3.vertex_values)
499	Q = quantity1quantity2 + 3quantity3
500	assert allclose(Q.vertex_values,
501	quantity1.vertex_values * quantity2.vertex_values +
502	3*quantity3.vertex_values)
503	Q = quantity1quantity2 + 3quantity3 + 5.0
504	assert allclose(Q.vertex_values,
505	quantity1.vertex_values * quantity2.vertex_values +
506	3*quantity3.vertex_values + 5)
507
508	Q = quantity1*quantity2 - quantity3
509	assert allclose(Q.vertex_values,
510	quantity1.vertex_values * quantity2.vertex_values -
511	quantity3.vertex_values)
512	Q = 1.5quantity1quantity2 - 3*quantity3 + 5.0
513	assert allclose(Q.vertex_values,
514	1.5quantity1.vertex_values quantity2.vertex_values -
515	3*quantity3.vertex_values + 5)
516
517	#Try combining quantities and arrays and scalars
518	Q = 1.5quantity1quantity2.vertex_values -\
519	3*quantity3.vertex_values + 5.0
520	assert allclose(Q.vertex_values,
521	1.5quantity1.vertex_values quantity2.vertex_values -
522	3*quantity3.vertex_values + 5)
523
524
525
526	def test_gradient(self):
527	quantity = Conserved_quantity(self.mesh4)
528
529	#Set up for a gradient of (2,0) at mid triangle
530	quantity.set_values([2.0, 4.0, 6.0, 2.0],
531	location = 'centroids')
532
533
534	#Gradients
535	a, b = quantity.compute_gradients()
536
537	#print self.mesh4.centroid_coordinates
538	#print a, b
539
540	#The central triangle (1)
541	#(using standard gradient based on neigbours controid values)
542	assert allclose(a[1], 2.0)
543	assert allclose(b[1], 0.0)
544
545
546	#Left triangle (0) using two point gradient
547	#q0 = q1 + a(x0-x1) + b(y0-y1) <=>
548	#2 = 4 + a(-2/3) + b(-2/3)
549	assert allclose(a[0] + b[0], 3)
550	#From orthogonality (a(y0-y1) + b(x0-x1) == 0)
551	assert allclose(a[0] - b[0], 0)
552
553
554	#Right triangle (2) using two point gradient
555	#q2 = q1 + a(x2-x1) + b(y2-y1) <=>
556	#6 = 4 + a(4/3) + b(-2/3)
557	assert allclose(2*a[2] - b[2], 3)
558	#From orthogonality (a(y1-y2) + b(x2-x1) == 0)
559	assert allclose(a[2] + 2*b[2], 0)
560
561
562	#Top triangle (3) using two point gradient
563	#q3 = q1 + a(x3-x1) + b(y3-y1) <=>
564	#2 = 4 + a(-2/3) + b(4/3)
565	assert allclose(a[3] - 2*b[3], 3)
566	#From orthogonality (a(y1-y3) + b(x3-x1) == 0)
567	assert allclose(2*a[3] + b[3], 0)
568
569
570
571	#print a, b
572	quantity.extrapolate_second_order()
573
574	#Apply q(x,y) = qc + a(x-xc) + b(y-yc)
575	assert allclose(quantity.vertex_values[0,:], [3., 0., 3.])
576	assert allclose(quantity.vertex_values[1,:], [4./3, 16./3, 16./3])
577
578
579	#a = 1.2, b=-0.6
580	#q(4,0) = 6 + a(4 - 8/3) + b(-2/3)
581	assert allclose(quantity.vertex_values[2,2], 8)
582
583
584	def test_second_order_extrapolation2(self):
585	quantity = Conserved_quantity(self.mesh4)
586
587	#Set up for a gradient of (3,1), f(x) = 3x+y
588	quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
589	location = 'centroids')
590
591	#Gradients
592	a, b = quantity.compute_gradients()
593
594	#print a, b
595
596	assert allclose(a[1], 3.0)
597	assert allclose(b[1], 1.0)
598
599	#Work out the others
600
601	quantity.extrapolate_second_order()
602
603	#print quantity.vertex_values
604	assert allclose(quantity.vertex_values[1,0], 2.0)
605	assert allclose(quantity.vertex_values[1,1], 6.0)
606	assert allclose(quantity.vertex_values[1,2], 8.0)
607
608
609
610	def test_first_order_extrapolator(self):
611	quantity = Conserved_quantity(self.mesh4)
612
613	#Test centroids
614	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
615	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
616
617	#Extrapolate
618	quantity.extrapolate_first_order()
619
620	#Check vertices but not edge values
621	assert allclose(quantity.vertex_values,
622	[[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
623
624
625	def test_second_order_extrapolator(self):
626	quantity = Conserved_quantity(self.mesh4)
627
628	#Set up for a gradient of (3,0) at mid triangle
629	quantity.set_values([2.0, 4.0, 8.0, 2.0],
630	location = 'centroids')
631
632
633
634	quantity.extrapolate_second_order()
635	quantity.limit()
636
637
638	#Assert that central triangle is limited by neighbours
639	assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
640	assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
641
642	assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
643	assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
644
645	assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
646	assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
647
648
649	#Assert that quantities are conserved
650	from Numeric import sum
651	for k in range(quantity.centroid_values.shape[0]):
652	assert allclose (quantity.centroid_values[k],
653	sum(quantity.vertex_values[k,:])/3)
654
655
656
657
658
659	def test_limiter(self):
660	quantity = Conserved_quantity(self.mesh4)
661
662	#Create a deliberate overshoot (e.g. from gradient computation)
663	quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
664
665
666	#Limit
667	quantity.limit()
668
669	#Assert that central triangle is limited by neighbours
670	assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
671	assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
672
673	assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
674	assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
675
676	assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
677	assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
678
679
680
681	#Assert that quantities are conserved
682	from Numeric import sum
683	for k in range(quantity.centroid_values.shape[0]):
684	assert allclose (quantity.centroid_values[k],
685	sum(quantity.vertex_values[k,:])/3)
686
687
688	def test_limiter2(self):
689	"""Taken from test_shallow_water
690	"""
691	quantity = Conserved_quantity(self.mesh4)
692
693	#Test centroids
694	quantity.set_values([2.,4.,8.,2.], location = 'centroids')
695	assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
696
697
698	#Extrapolate
699	quantity.extrapolate_second_order()
700
701	assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
702
703	#Limit
704	quantity.limit()
705
706
707	assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
708
709
710	#Assert that quantities are conserved
711	from Numeric import sum
712	for k in range(quantity.centroid_values.shape[0]):
713	assert allclose (quantity.centroid_values[k],
714	sum(quantity.vertex_values[k,:])/3)
715
716
717
718
719
720	def test_distribute_first_order(self):
721	quantity = Conserved_quantity(self.mesh4)
722
723	#Test centroids
724	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
725	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
726
727
728	#Extrapolate
729	quantity.extrapolate_first_order()
730
731	#Interpolate
732	quantity.interpolate_from_vertices_to_edges()
733
734	assert allclose(quantity.vertex_values,
735	[[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
736	assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
737	[3,3,3], [4, 4, 4]])
738
739
740	def test_distribute_second_order(self):
741	quantity = Conserved_quantity(self.mesh4)
742
743	#Test centroids
744	quantity.set_values([2.,4.,8.,2.], location = 'centroids')
745	assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
746
747
748	#Extrapolate
749	quantity.extrapolate_second_order()
750
751	assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
752
753
754	def test_update_explicit(self):
755	quantity = Conserved_quantity(self.mesh4)
756
757	#Test centroids
758	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
759	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
760
761	#Set explicit_update
762	quantity.explicit_update = array( [1.,1.,1.,1.] )
763
764	#Update with given timestep
765	quantity.update(0.1)
766
767	x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
768	assert allclose( quantity.centroid_values, x)
769
770	def test_update_semi_implicit(self):
771	quantity = Conserved_quantity(self.mesh4)
772
773	#Test centroids
774	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
775	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
776
777	#Set semi implicit update
778	quantity.semi_implicit_update = array([1.,1.,1.,1.])
779
780	#Update with given timestep
781	timestep = 0.1
782	quantity.update(timestep)
783
784	sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
785	denom = ones(4, Float)-timestep*sem
786
787	x = array([1, 2, 3, 4])/denom
788	assert allclose( quantity.centroid_values, x)
789
790
791	def test_both_updates(self):
792	quantity = Conserved_quantity(self.mesh4)
793
794	#Test centroids
795	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
796	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
797
798	#Set explicit_update
799	quantity.explicit_update = array( [4.,3.,2.,1.] )
800
801	#Set semi implicit update
802	quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
803
804	#Update with given timestep
805	timestep = 0.1
806	quantity.update(0.1)
807
808	sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
809	denom = ones(4, Float)-timestep*sem
810
811	x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
812	x /= denom
813	assert allclose( quantity.centroid_values, x)
814
815
816
817
818	#Test smoothing
819	def test_smoothing(self):
820
821	from mesh_factory import rectangular
822	from shallow_water import Domain, Transmissive_boundary
823	from Numeric import zeros, Float
824	from util import mean
825
826	#Create basic mesh
827	points, vertices, boundary = rectangular(2, 2)
828
829	#Create shallow water domain
830	domain = Domain(points, vertices, boundary)
831	domain.default_order=2
832	domain.reduction = mean
833
834
835	#Set some field values
836	domain.set_quantity('elevation', lambda x,y: x)
837	domain.set_quantity('friction', 0.03)
838
839
840	######################
841	# Boundary conditions
842	B = Transmissive_boundary(domain)
843	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
844
845
846	######################
847	#Initial condition - with jumps
848
849	bed = domain.quantities['elevation'].vertex_values
850	stage = zeros(bed.shape, Float)
851
852	h = 0.03
853	for i in range(stage.shape[0]):
854	if i % 2 == 0:
855	stage[i,:] = bed[i,:] + h
856	else:
857	stage[i,:] = bed[i,:]
858
859	domain.set_quantity('stage', stage)
860
861	stage = domain.quantities['stage']
862
863	#Get smoothed stage
864	A, V = stage.get_vertex_values(xy=False, smooth=True)
865	Q = stage.vertex_values
866
867
868	assert A.shape[0] == 9
869	assert V.shape[0] == 8
870	assert V.shape[1] == 3
871
872	#First four points
873	assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
874	assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
875	assert allclose(A[2], Q[3,0])
876	assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
877
878	#Center point
879	assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
880	Q[5,0] + Q[6,2] + Q[7,1])/6)
881
882
883	#Check V
884	assert allclose(V[0,:], [3,4,0])
885	assert allclose(V[1,:], [1,0,4])
886	assert allclose(V[2,:], [4,5,1])
887	assert allclose(V[3,:], [2,1,5])
888	assert allclose(V[4,:], [6,7,3])
889	assert allclose(V[5,:], [4,3,7])
890	assert allclose(V[6,:], [7,8,4])
891	assert allclose(V[7,:], [5,4,8])
892
893	#Get smoothed stage with XY
894	X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
895
896	assert allclose(A, A1)
897	assert allclose(V, V1)
898
899	#Check XY
900	assert allclose(X[4], 0.5)
901	assert allclose(Y[4], 0.5)
902
903	assert allclose(X[7], 1.0)
904	assert allclose(Y[7], 0.5)
905
906
907
908
909	def test_vertex_values_no_smoothing(self):
910
911	from mesh_factory import rectangular
912	from shallow_water import Domain, Transmissive_boundary
913	from Numeric import zeros, Float
914	from util import mean
915
916
917	#Create basic mesh
918	points, vertices, boundary = rectangular(2, 2)
919
920	#Create shallow water domain
921	domain = Domain(points, vertices, boundary)
922	domain.default_order=2
923	domain.reduction = mean
924
925
926	#Set some field values
927	domain.set_quantity('elevation', lambda x,y: x)
928	domain.set_quantity('friction', 0.03)
929
930
931	######################
932	#Initial condition - with jumps
933
934	bed = domain.quantities['elevation'].vertex_values
935	stage = zeros(bed.shape, Float)
936
937	h = 0.03
938	for i in range(stage.shape[0]):
939	if i % 2 == 0:
940	stage[i,:] = bed[i,:] + h
941	else:
942	stage[i,:] = bed[i,:]
943
944	domain.set_quantity('stage', stage)
945
946	#Get stage
947	stage = domain.quantities['stage']
948	A, V = stage.get_vertex_values(xy=False, smooth=False)
949	Q = stage.vertex_values.flat
950
951	for k in range(8):
952	assert allclose(A[k], Q[k])
953
954
955	for k in range(8):
956	assert V[k, 0] == 3*k
957	assert V[k, 1] == 3*k+1
958	assert V[k, 2] == 3*k+2
959
960
961
962	X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
963
964
965	assert allclose(A, A1)
966	assert allclose(V, V1)
967
968	#Check XY
969	assert allclose(X[1], 0.5)
970	assert allclose(Y[1], 0.5)
971	assert allclose(X[4], 0.0)
972	assert allclose(Y[4], 0.0)
973	assert allclose(X[12], 1.0)
974	assert allclose(Y[12], 0.0)
975
976
977
978	def set_array_values_by_index(self):
979
980	from mesh_factory import rectangular
981	from shallow_water import Domain
982	from Numeric import zeros, Float
983
984	#Create basic mesh
985	points, vertices, boundary = rectangular(1, 1)
986
987	#Create shallow water domain
988	domain = Domain(points, vertices, boundary)
989	#print "domain.number_of_elements ",domain.number_of_elements
990	quantity = Quantity(domain,[[1,1,1],[2,2,2]])
991	value = [7]
992	indices = [1]
993	quantity.set_array_values_by_index(value,
994	location = 'centroids',
995	indices = indices)
996	#print "quantity.centroid_values",quantity.centroid_values
997
998	assert allclose(quantity.centroid_values, [1,7])
999
1000	quantity.set_array_values([15,20,25], indices = indices)
1001	assert allclose(quantity.centroid_values, [1,20])
1002
1003	quantity.set_array_values([15,20,25], indices = indices)
1004	assert allclose(quantity.centroid_values, [1,20])
1005
1006	def test_setting_some_vertex_values(self):
1007	"""
1008	set values based on triangle lists.
1009	"""
1010	from mesh_factory import rectangular
1011	from shallow_water import Domain
1012	from Numeric import zeros, Float
1013
1014	#Create basic mesh
1015	points, vertices, boundary = rectangular(1, 3)
1016	#print "vertices",vertices
1017	#Create shallow water domain
1018	domain = Domain(points, vertices, boundary)
1019	#print "domain.number_of_elements ",domain.number_of_elements
1020	quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1021	[4,4,4],[5,5,5],[6,6,6]])
1022	value = [7]
1023	indices = [1]
1024	quantity.set_values(value,
1025	location = 'centroids',
1026	indices = indices)
1027	#print "quantity.centroid_values",quantity.centroid_values
1028	assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1029
1030	value = [[15,20,25]]
1031	quantity.set_values(value, indices = indices)
1032	#print "1 quantity.vertex_values",quantity.vertex_values
1033	assert allclose(quantity.vertex_values[1], value[0])
1034
1035
1036	#print "quantity",quantity.vertex_values
1037	values = [10,100,50]
1038	quantity.set_values(values, indices = [0,1,5], location = 'centroids')
1039	#print "2 quantity.vertex_values",quantity.vertex_values
1040	assert allclose(quantity.vertex_values[0], [10,10,10])
1041	assert allclose(quantity.vertex_values[5], [50,50,50])
1042	#quantity.interpolate()
1043	#print "quantity.centroid_values",quantity.centroid_values
1044	assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
1045
1046
1047	quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1048	[4,4,4],[5,5,5],[6,6,6]])
1049	values = [10,100,50]
1050	#this will be per unique vertex, indexing the vertices
1051	#print "quantity.vertex_values",quantity.vertex_values
1052	quantity.set_values(values, indices = [0,1,5])
1053	#print "quantity.vertex_values",quantity.vertex_values
1054	assert allclose(quantity.vertex_values[0], [1,50,10])
1055	assert allclose(quantity.vertex_values[5], [6,6,6])
1056	assert allclose(quantity.vertex_values[1], [100,10,50])
1057
1058	quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1059	[4,4,4],[5,5,5],[6,6,6]])
1060	values = [[31,30,29],[400,400,400],[1000,999,998]]
1061	quantity.set_values(values, indices = [3,3,5])
1062	quantity.interpolate()
1063	assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
1064
1065	values = [[1,1,1],[2,2,2],[3,3,3],
1066	[4,4,4],[5,5,5],[6,6,6]]
1067	quantity.set_values(values)
1068
1069	# testing the standard set values by vertex
1070	# indexed by vertex_id in general_mesh.coordinates
1071	values = [0,1,2,3,4,5,6,7]
1072
1073	quantity.set_values(values)
1074	#print "1 quantity.vertex_values",quantity.vertex_values
1075	assert allclose(quantity.vertex_values,[[ 4., 5., 0.],
1076	[ 1., 0., 5.],
1077	[ 5., 6., 1.],
1078	[ 2., 1., 6.],
1079	[ 6., 7., 2.],
1080	[ 3., 2., 7.]])
1081
1082	def test_setting_unique_vertex_values(self):
1083	"""
1084	set values based on unique_vertex lists.
1085	"""
1086	from mesh_factory import rectangular
1087	from shallow_water import Domain
1088	from Numeric import zeros, Float
1089
1090	#Create basic mesh
1091	points, vertices, boundary = rectangular(1, 3)
1092	#print "vertices",vertices
1093	#Create shallow water domain
1094	domain = Domain(points, vertices, boundary)
1095	#print "domain.number_of_elements ",domain.number_of_elements
1096	quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1097	[4,4,4],[5,5,5]])
1098	value = 7
1099	indices = [1,5]
1100	quantity.set_values(value,
1101	location = 'unique vertices',
1102	indices = indices)
1103	#print "quantity.centroid_values",quantity.centroid_values
1104	assert allclose(quantity.vertex_values[0], [0,7,0])
1105	assert allclose(quantity.vertex_values[1], [7,1,7])
1106	assert allclose(quantity.vertex_values[2], [7,2,7])
1107
1108
1109	def test_get_values(self):
1110	"""
1111	get values based on triangle lists.
1112	"""
1113	from mesh_factory import rectangular
1114	from shallow_water import Domain
1115	from Numeric import zeros, Float
1116
1117	#Create basic mesh
1118	points, vertices, boundary = rectangular(1, 3)
1119
1120	#print "points",points
1121	#print "vertices",vertices
1122	#print "boundary",boundary
1123
1124	#Create shallow water domain
1125	domain = Domain(points, vertices, boundary)
1126	#print "domain.number_of_elements ",domain.number_of_elements
1127	quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1128	[4,4,4],[5,5,5]])
1129
1130	#print "quantity.get_values(location = 'unique vertices')", \
1131	# quantity.get_values(location = 'unique vertices')
1132
1133	#print "quantity.get_values(location = 'unique vertices')", \
1134	# quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
1135	# location = 'unique vertices')
1136
1137	answer = [0.5,2,4,5,0,1,3,4.5]
1138	assert allclose(answer,
1139	quantity.get_values(location = 'unique vertices'))
1140
1141	indices = [0,5,3]
1142	answer = [0.5,1,5]
1143	assert allclose(answer,
1144	quantity.get_values(indices=indices, \
1145	location = 'unique vertices'))
1146	#print "quantity.centroid_values",quantity.centroid_values
1147	#print "quantity.get_values(location = 'centroids') ",\
1148	# quantity.get_values(location = 'centroids')
1149
1150	def test_getting_some_vertex_values(self):
1151	"""
1152	get values based on triangle lists.
1153	"""
1154	from mesh_factory import rectangular
1155	from shallow_water import Domain
1156	from Numeric import zeros, Float
1157
1158	#Create basic mesh
1159	points, vertices, boundary = rectangular(1, 3)
1160
1161	#print "points",points
1162	#print "vertices",vertices
1163	#print "boundary",boundary
1164
1165	#Create shallow water domain
1166	domain = Domain(points, vertices, boundary)
1167	#print "domain.number_of_elements ",domain.number_of_elements
1168	quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1169	[4,4,4],[5,5,5],[6,6,6]])
1170	value = [7]
1171	indices = [1]
1172	quantity.set_values(value,
1173	location = 'centroids',
1174	indices = indices)
1175	#print "quantity.centroid_values",quantity.centroid_values
1176	#print "quantity.get_values(location = 'centroids') ",\
1177	# quantity.get_values(location = 'centroids')
1178	assert allclose(quantity.centroid_values,
1179	quantity.get_values(location = 'centroids'))
1180
1181
1182	value = [[15,20,25]]
1183	quantity.set_values(value, indices = indices)
1184	#print "1 quantity.vertex_values",quantity.vertex_values
1185	assert allclose(quantity.vertex_values, quantity.get_values())
1186
1187	assert allclose(quantity.edge_values,
1188	quantity.get_values(location = 'edges'))
1189
1190	# get a subset of elements
1191	subset = quantity.get_values(location='centroids', indices=[0,5])
1192	answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
1193	assert allclose(subset, answer)
1194
1195
1196	subset = quantity.get_values(location='edges', indices=[0,5])
1197	answer = [quantity.edge_values[0],quantity.edge_values[5]]
1198	#print "subset",subset
1199	#print "answer",answer
1200	assert allclose(subset, answer)
1201
1202	subset = quantity.get_values( indices=[1,5])
1203	answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
1204	#print "subset",subset
1205	#print "answer",answer
1206	assert allclose(subset, answer)
1207
1208
1209
1210	#-------------------------------------------------------------
1211	if __name__ == "__main__":
1212	suite = unittest.makeSuite(Test_Quantity,'test')
1213	#print "restricted test"
1214	#suite = unittest.makeSuite(Test_Quantity,'test_set_values_func')
1215	runner = unittest.TextTestRunner()
1216	runner.run(suite)

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

Download in other formats: