Context Navigation

source: inundation/ga/storm_surge/pyvolution/test_shallow_water.py @ 668

Last change on this file since 668 was 648, checked in by ole, 20 years ago
Refactoring boundary structure to allow None Boundary object for internal boundaries.
File size: 81.7 KB

Rev	Line
[218]	1	#!/usr/bin/env python
	2
	3	import unittest
	4	from math import sqrt, pi
	5
	6	from config import g, epsilon
	7	from Numeric import allclose, array, zeros, ones, Float
	8	from shallow_water import *
	9
[609]	10
	11
	12	#Variable windfield implemented using functions
[623]	13	def speed(t,x,y):
[609]	14	"""Large speeds halfway between center and edges
	15	Low speeds at center and edges
	16	"""
	17
	18	from math import sqrt, exp, cos, pi
[610]	19
	20	N = len(x)
	21	s = 0*x #New array
[609]	22
[610]	23	for k in range(N):
	24
	25	r = sqrt(x[k]2 + y[k]2)
[609]	26
[610]	27	factor = exp( -(r-0.15)**2 )
[609]	28
[610]	29	s[k] = 4000 * factor * (cos(t2pi/150) + 2)
[609]	30
[610]	31	return s
[609]	32
[610]	33
[623]	34	def scalar_func(t,x,y):
[612]	35	"""Function that returns a scalar.
	36	Used to test error message when Numeric array is expected
	37	"""
	38
	39	return 17.7
	40
	41
[623]	42	def angle(t,x,y):
[609]	43	"""Rotating field
	44	"""
	45	from math import sqrt, atan, pi
	46
	47
[610]	48	N = len(x)
	49	a = 0*x #New array
[609]	50
[610]	51	for k in range(N):
	52	r = sqrt(x[k]2 + y[k]2)
[609]	53
[610]	54	angle = atan(y[k]/x[k])
[609]	55
[610]	56	if x[k] < 0:
	57	angle+=pi #atan in ]-pi/2; pi/2[
[609]	58
[610]	59	#Take normal direction
	60	angle -= pi/2
	61
	62	#Ensure positive radians
	63	if angle < 0:
	64	angle += 2*pi
[609]	65
[610]	66	a[k] = angle/pi*180
[218]	67
[610]	68	return a
	69
	70
[218]	71	class TestCase(unittest.TestCase):
	72	def setUp(self):
	73	pass
	74
	75	def tearDown(self):
	76	pass
	77
[240]	78	def test_rotate(self):
	79	normal = array([0.0,-1.0])
[218]	80
[240]	81	q = array([1.0,2.0,3.0])
	82
	83	r = rotate(q, normal, direction = 1)
	84	assert r[0] == 1
	85	assert r[1] == -3
	86	assert r[2] == 2
	87
	88	w = rotate(r, normal, direction = -1)
	89	assert allclose(w, q)
	90
	91
	92
[218]	93	def test_flux_zero_case(self):
	94	ql = zeros( 3, Float )
	95	qr = zeros( 3, Float )
	96	normal = zeros( 2, Float )
	97	zl = zr = 0.
	98	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	99
	100	assert allclose(flux, [0,0,0])
	101	assert max_speed == 0.
	102
	103	def test_flux_constants(self):
	104	w = 2.0
	105
	106	normal = array([1.,0])
	107	ql = array([w, 0, 0])
	108	qr = array([w, 0, 0])
	109	zl = zr = 0.
	110	h = w - (zl+zr)/2
	111
	112	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	113
	114	assert allclose(flux, [0., 0.5gh**2, 0.])
	115	assert max_speed == sqrt(g*h)
	116
	117
	118	def test_flux1(self):
	119	#Use data from previous version of pyvolution
	120	normal = array([1.,0])
	121	ql = array([-0.2, 2, 3])
	122	qr = array([-0.2, 2, 3])
	123	zl = zr = -0.5
	124	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	125
	126	assert allclose(flux, [2.,13.77433333, 20.])
	127	assert allclose(max_speed, 8.38130948661)
	128
	129
	130	def test_flux2(self):
	131	#Use data from previous version of pyvolution
	132	normal = array([0., -1.])
	133	ql = array([-0.075, 2, 3])
	134	qr = array([-0.075, 2, 3])
	135	zl = zr = -0.375
	136	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	137
	138	assert allclose(flux, [-3.,-20.0, -30.441])
	139	assert allclose(max_speed, 11.7146428199)
	140
	141	def test_flux3(self):
	142	#Use data from previous version of pyvolution
	143	normal = array([-sqrt(2)/2, sqrt(2)/2])
	144	ql = array([-0.075, 2, 3])
	145	qr = array([-0.075, 2, 3])
	146	zl = zr = -0.375
	147	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	148
	149	assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
	150	assert allclose(max_speed, 4.0716654239)
	151
	152	def test_flux4(self):
	153	#Use data from previous version of pyvolution
	154	normal = array([-sqrt(2)/2, sqrt(2)/2])
	155	ql = array([-0.34319278, 0.10254161, 0.07273855])
	156	qr = array([-0.30683287, 0.1071986, 0.05930515])
	157	zl = zr = -0.375
	158	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	159
	160	assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
	161	assert allclose(max_speed, 1.31414103233)
	162
	163	def test_sw_domain_simple(self):
	164	a = [0.0, 0.0]
	165	b = [0.0, 2.0]
	166	c = [2.0,0.0]
	167	d = [0.0, 4.0]
	168	e = [2.0, 2.0]
	169	f = [4.0,0.0]
	170
	171	points = [a, b, c, d, e, f]
	172	#bac, bce, ecf, dbe, daf, dae
	173	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
	174
	175	domain = Domain(points, vertices)
	176	domain.check_integrity()
	177
	178	for name in ['level', 'xmomentum', 'ymomentum',
	179	'elevation', 'friction']:
	180	assert domain.quantities.has_key(name)
	181
	182
	183	assert domain.get_conserved_quantities(0, edge=1) == 0.
	184
	185
	186	def test_boundary_conditions(self):
	187
	188	a = [0.0, 0.0]
	189	b = [0.0, 2.0]
	190	c = [2.0,0.0]
	191	d = [0.0, 4.0]
	192	e = [2.0, 2.0]
	193	f = [4.0,0.0]
	194
	195	points = [a, b, c, d, e, f]
	196	#bac, bce, ecf, dbe
	197	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
	198	boundary = { (0, 0): 'Third',
	199	(0, 2): 'First',
	200	(2, 0): 'Second',
	201	(2, 1): 'Second',
	202	(3, 1): 'Second',
	203	(3, 2): 'Third'}
	204
	205
	206	domain = Domain(points, vertices, boundary)
	207	domain.check_integrity()
	208
	209
	210	domain.set_quantity('level', [[1,2,3], [5,5,5],
	211	[0,0,9], [-6, 3, 3]])
	212
	213	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
	214	[3,3,3], [4, 4, 4]])
	215
	216	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
	217	[30,30,30], [40, 40, 40]])
	218
	219
	220	D = Dirichlet_boundary([5,2,1])
	221	T = Transmissive_boundary(domain)
	222	R = Reflective_boundary(domain)
	223	domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
	224
	225	domain.update_boundary()
	226
	227	#Level
	228	assert domain.quantities['level'].boundary_values[0] == 2.5
	229	assert domain.quantities['level'].boundary_values[0] ==\
	230	domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
	231	assert domain.quantities['level'].boundary_values[1] == 5. #Dirichlet
	232	assert domain.quantities['level'].boundary_values[2] ==\
	233	domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
	234	assert domain.quantities['level'].boundary_values[3] ==\
	235	domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
	236	assert domain.quantities['level'].boundary_values[4] ==\
	237	domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
	238	assert domain.quantities['level'].boundary_values[5] ==\
	239	domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
	240
	241	#Xmomentum
	242	assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
	243	assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
	244	assert domain.quantities['xmomentum'].boundary_values[2] ==\
	245	domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
	246	assert domain.quantities['xmomentum'].boundary_values[3] ==\
	247	domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
	248	assert domain.quantities['xmomentum'].boundary_values[4] ==\
	249	domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
	250	assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 #Reflective
	251
	252	#Ymomentum
	253	assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
	254	assert domain.quantities['ymomentum'].boundary_values[1] == 1. #Dirichlet
	255	assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
	256	assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
	257	assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
	258	assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
	259
	260
[646]	261	def test_boundary_conditionsII(self):
[218]	262
[646]	263	a = [0.0, 0.0]
	264	b = [0.0, 2.0]
	265	c = [2.0,0.0]
	266	d = [0.0, 4.0]
	267	e = [2.0, 2.0]
	268	f = [4.0,0.0]
[218]	269
[646]	270	points = [a, b, c, d, e, f]
	271	#bac, bce, ecf, dbe
	272	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
	273	boundary = { (0, 0): 'Third',
	274	(0, 2): 'First',
	275	(2, 0): 'Second',
	276	(2, 1): 'Second',
	277	(3, 1): 'Second',
	278	(3, 2): 'Third',
	279	(0, 1): 'Internal'}
	280
	281
	282	domain = Domain(points, vertices, boundary)
	283	domain.check_integrity()
	284
	285
	286	domain.set_quantity('level', [[1,2,3], [5,5,5],
	287	[0,0,9], [-6, 3, 3]])
	288
	289	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
	290	[3,3,3], [4, 4, 4]])
	291
	292	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
	293	[30,30,30], [40, 40, 40]])
	294
	295
	296	D = Dirichlet_boundary([5,2,1])
	297	T = Transmissive_boundary(domain)
	298	R = Reflective_boundary(domain)
[648]	299	domain.set_boundary( {'First': D, 'Second': T,
	300	'Third': R, 'Internal': None})
	301
[646]	302	domain.update_boundary()
[648]	303	domain.check_integrity()
[646]	304
	305
[240]	306	def test_compute_fluxes0(self):
[218]	307	#Do a full triangle and check that fluxes cancel out for
	308	#the constant level case
	309
	310	a = [0.0, 0.0]
	311	b = [0.0, 2.0]
	312	c = [2.0,0.0]
	313	d = [0.0, 4.0]
	314	e = [2.0, 2.0]
	315	f = [4.0,0.0]
	316
	317	points = [a, b, c, d, e, f]
	318	#bac, bce, ecf, dbe
	319	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	320
	321	domain = Domain(points, vertices)
	322	domain.set_quantity('level', [[2,2,2], [2,2,2],
	323	[2,2,2], [2,2,2]])
	324	domain.check_integrity()
	325
	326	zl=zr=0. #Assume flat bed
	327
	328	#Flux across right edge of volume 1
	329	normal = domain.get_normal(1,0)
	330	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
	331	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
	332	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
	333
	334	#Check that flux seen from other triangles is inverse
	335	tmp = qr; qr=ql; ql=tmp
	336	normal = domain.get_normal(2,2)
	337	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	338	assert allclose(flux + flux0, 0.)
	339
	340	#Flux across upper edge of volume 1
	341	normal = domain.get_normal(1,1)
	342	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
	343	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
	344	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
	345
	346	#Check that flux seen from other triangles is inverse
	347	tmp = qr; qr=ql; ql=tmp
	348	normal = domain.get_normal(3,0)
	349	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	350	assert allclose(flux + flux1, 0.)
	351
	352	#Flux across lower left hypotenuse of volume 1
	353	normal = domain.get_normal(1,2)
	354	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
	355	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
	356	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
	357
	358	#Check that flux seen from other triangles is inverse
	359	tmp = qr; qr=ql; ql=tmp
	360	normal = domain.get_normal(0,1)
	361	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	362	assert allclose(flux + flux2, 0.)
	363
	364
	365	#Scale by edgelengths, add up anc check that total flux is zero
	366	e0 = domain.edgelengths[1, 0]
	367	e1 = domain.edgelengths[1, 1]
	368	e2 = domain.edgelengths[1, 2]
	369
	370	assert allclose(e0flux0+e1flux1+e2*flux2, 0.)
	371
	372	#Now check that compute_flux yields zeros as well
	373	domain.compute_fluxes()
	374
	375	for name in ['level', 'xmomentum', 'ymomentum']:
	376	#print name, domain.quantities[name].explicit_update
	377	assert allclose(domain.quantities[name].explicit_update[1], 0)
	378
	379
	380
	381	def test_compute_fluxes1(self):
	382	#Use values from previous version
	383
	384	a = [0.0, 0.0]
	385	b = [0.0, 2.0]
	386	c = [2.0,0.0]
	387	d = [0.0, 4.0]
	388	e = [2.0, 2.0]
	389	f = [4.0,0.0]
	390
	391	points = [a, b, c, d, e, f]
	392	#bac, bce, ecf, dbe
	393	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	394
	395	domain = Domain(points, vertices)
	396	val0 = 2.+2.0/3
	397	val1 = 4.+4.0/3
	398	val2 = 8.+2.0/3
	399	val3 = 2.+8.0/3
	400
	401	domain.set_quantity('level', [[val0, val0, val0], [val1, val1, val1],
	402	[val2, val2, val2], [val3, val3, val3]])
	403	domain.check_integrity()
	404
	405	zl=zr=0. #Assume flat bed
	406
	407	#Flux across right edge of volume 1
	408	normal = domain.get_normal(1,0)
	409	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
	410	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
	411	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
	412	assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
	413	assert allclose(max_speed, 9.21592824046)
	414
	415
	416	#Flux across upper edge of volume 1
	417	normal = domain.get_normal(1,1)
	418	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
	419	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
	420	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
	421	assert allclose(flux1, [2.4098563, 0., 123.04444444])
	422	assert allclose(max_speed, 7.22956891292)
	423
	424	#Flux across lower left hypotenuse of volume 1
	425	normal = domain.get_normal(1,2)
	426	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
	427	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
	428	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
	429
	430	assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
	431	assert allclose(max_speed, 7.22956891292)
	432
	433	#Scale, add up and check that compute_fluxes is correct for vol 1
	434	e0 = domain.edgelengths[1, 0]
	435	e1 = domain.edgelengths[1, 1]
	436	e2 = domain.edgelengths[1, 2]
	437
	438	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
	439	assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
	440
	441
	442	domain.compute_fluxes()
	443
	444	#assert allclose(total_flux, domain.explicit_update[1,:])
	445	for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
	446	assert allclose(total_flux[i],
	447	domain.quantities[name].explicit_update[1])
	448
	449	#assert allclose(domain.explicit_update, [
	450	# [0., -69.68888889, -69.68888889],
	451	# [-0.68218178, -166.6, -35.93333333],
	452	# [-111.77316251, 69.68888889, 0.],
	453	# [-35.68522449, 0., 69.68888889]])
	454
	455	assert allclose(domain.quantities['level'].explicit_update,
	456	[0., -0.68218178, -111.77316251, -35.68522449])
	457	assert allclose(domain.quantities['xmomentum'].explicit_update,
	458	[-69.68888889, -166.6, 69.68888889, 0])
	459	assert allclose(domain.quantities['ymomentum'].explicit_update,
	460	[-69.68888889, -35.93333333, 0., 69.68888889])
	461
	462
	463	#assert allclose(domain.quantities[name].explicit_update
	464
	465
	466
	467
	468
	469	def test_compute_fluxes2(self):
	470	#Random values, incl momentum
	471
	472	a = [0.0, 0.0]
	473	b = [0.0, 2.0]
	474	c = [2.0,0.0]
	475	d = [0.0, 4.0]
	476	e = [2.0, 2.0]
	477	f = [4.0,0.0]
	478
	479	points = [a, b, c, d, e, f]
	480	#bac, bce, ecf, dbe
	481	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	482
	483	domain = Domain(points, vertices)
	484	val0 = 2.+2.0/3
	485	val1 = 4.+4.0/3
	486	val2 = 8.+2.0/3
	487	val3 = 2.+8.0/3
	488
	489	zl=zr=0 #Assume flat zero bed
	490
	491	domain.set_quantity('elevation', zl*ones( (4,3) ))
	492
	493
	494	domain.set_quantity('level', [[val0, val0-1, val0-2],
	495	[val1, val1+1, val1],
	496	[val2, val2-2, val2],
	497	[val3-0.5, val3, val3]])
	498
	499	domain.set_quantity('xmomentum',
	500	[[1, 2, 3], [3, 4, 5],
	501	[1, -1, 0], [0, -2, 2]])
	502
	503	domain.set_quantity('ymomentum',
	504	[[1, -1, 0], [0, -3, 2],
	505	[0, 1, 0], [-1, 2, 2]])
	506
	507
	508	domain.check_integrity()
	509
	510
	511
	512	#Flux across right edge of volume 1
	513	normal = domain.get_normal(1,0)
	514	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
	515	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
	516	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
	517
	518	#Flux across upper edge of volume 1
	519	normal = domain.get_normal(1,1)
	520	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
	521	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
	522	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
	523
	524	#Flux across lower left hypotenuse of volume 1
	525	normal = domain.get_normal(1,2)
	526	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
	527	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
	528	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
	529
	530	#Scale, add up and check that compute_fluxes is correct for vol 1
	531	e0 = domain.edgelengths[1, 0]
	532	e1 = domain.edgelengths[1, 1]
	533	e2 = domain.edgelengths[1, 2]
	534
	535	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
	536
	537
	538	domain.compute_fluxes()
	539	for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
	540	assert allclose(total_flux[i],
	541	domain.quantities[name].explicit_update[1])
	542	#assert allclose(total_flux, domain.explicit_update[1,:])
	543
	544
	545	def test_compute_fluxes3(self):
	546	#Random values, incl momentum
	547
	548	a = [0.0, 0.0]
	549	b = [0.0, 2.0]
	550	c = [2.0,0.0]
	551	d = [0.0, 4.0]
	552	e = [2.0, 2.0]
	553	f = [4.0,0.0]
	554
	555	points = [a, b, c, d, e, f]
	556	#bac, bce, ecf, dbe
	557	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	558
	559	domain = Domain(points, vertices)
	560	val0 = 2.+2.0/3
	561	val1 = 4.+4.0/3
	562	val2 = 8.+2.0/3
	563	val3 = 2.+8.0/3
	564
	565	zl=zr=-3.75 #Assume constant bed (must be less than level)
	566	domain.set_quantity('elevation', zl*ones( (4,3) ))
	567
	568
	569	domain.set_quantity('level', [[val0, val0-1, val0-2],
	570	[val1, val1+1, val1],
	571	[val2, val2-2, val2],
	572	[val3-0.5, val3, val3]])
	573
	574	domain.set_quantity('xmomentum',
	575	[[1, 2, 3], [3, 4, 5],
	576	[1, -1, 0], [0, -2, 2]])
	577
	578	domain.set_quantity('ymomentum',
	579	[[1, -1, 0], [0, -3, 2],
	580	[0, 1, 0], [-1, 2, 2]])
	581
	582
	583	domain.check_integrity()
	584
	585
	586
	587	#Flux across right edge of volume 1
	588	normal = domain.get_normal(1,0)
	589	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
	590	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
	591	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
	592
	593	#Flux across upper edge of volume 1
	594	normal = domain.get_normal(1,1)
	595	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
	596	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
	597	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
	598
	599	#Flux across lower left hypotenuse of volume 1
	600	normal = domain.get_normal(1,2)
	601	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
	602	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
	603	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
	604
	605	#Scale, add up and check that compute_fluxes is correct for vol 1
	606	e0 = domain.edgelengths[1, 0]
	607	e1 = domain.edgelengths[1, 1]
	608	e2 = domain.edgelengths[1, 2]
	609
	610	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
	611
	612	domain.compute_fluxes()
	613	for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
	614	assert allclose(total_flux[i],
	615	domain.quantities[name].explicit_update[1])
	616
	617
	618
	619	def test_catching_negative_heights(self):
	620
	621	a = [0.0, 0.0]
	622	b = [0.0, 2.0]
	623	c = [2.0,0.0]
	624	d = [0.0, 4.0]
	625	e = [2.0, 2.0]
	626	f = [4.0,0.0]
	627
	628	points = [a, b, c, d, e, f]
	629	#bac, bce, ecf, dbe
	630	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	631
	632	domain = Domain(points, vertices)
	633	val0 = 2.+2.0/3
	634	val1 = 4.+4.0/3
	635	val2 = 8.+2.0/3
	636	val3 = 2.+8.0/3
	637
	638	zl=zr=4 #Too large
	639	domain.set_quantity('elevation', zl*ones( (4,3) ))
	640	domain.set_quantity('level', [[val0, val0-1, val0-2],
	641	[val1, val1+1, val1],
	642	[val2, val2-2, val2],
	643	[val3-0.5, val3, val3]])
	644
	645	#Should fail
	646	try:
	647	domain.check_integrity()
	648	except:
	649	pass
	650
	651	#####################################################
	652	def test_gravity(self):
	653	#Assuming no friction
	654
	655	from config import g
	656
	657	a = [0.0, 0.0]
	658	b = [0.0, 2.0]
	659	c = [2.0, 0.0]
	660	d = [0.0, 4.0]
	661	e = [2.0, 2.0]
	662	f = [4.0, 0.0]
	663
	664	points = [a, b, c, d, e, f]
	665	#bac, bce, ecf, dbe
	666	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	667
	668	domain = Domain(points, vertices)
	669
	670	#Set up for a gradient of (3,0) at mid triangle
	671	def slope(x, y):
	672	return 3*x
	673
	674	h = 0.1
	675	def level(x,y):
	676	return slope(x,y)+h
	677
	678	domain.set_quantity('elevation', slope)
	679	domain.set_quantity('level', level)
	680
	681	for name in domain.conserved_quantities:
	682	assert allclose(domain.quantities[name].explicit_update, 0)
	683	assert allclose(domain.quantities[name].semi_implicit_update, 0)
	684
	685	domain.compute_forcing_terms()
	686
	687	assert allclose(domain.quantities['level'].explicit_update, 0)
	688	assert allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
	689	assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
	690
	691
[609]	692	def test_manning_friction(self):
	693	from config import g
	694
	695	a = [0.0, 0.0]
	696	b = [0.0, 2.0]
	697	c = [2.0, 0.0]
	698	d = [0.0, 4.0]
	699	e = [2.0, 2.0]
	700	f = [4.0, 0.0]
[218]	701
[609]	702	points = [a, b, c, d, e, f]
	703	#bac, bce, ecf, dbe
	704	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	705
	706	domain = Domain(points, vertices)
	707
	708	#Set up for a gradient of (3,0) at mid triangle
	709	def slope(x, y):
	710	return 3*x
	711
	712	h = 0.1
	713	def level(x,y):
	714	return slope(x,y)+h
	715
	716	eta = 0.07
	717	domain.set_quantity('elevation', slope)
	718	domain.set_quantity('level', level)
	719	domain.set_quantity('friction', eta)
	720
	721	for name in domain.conserved_quantities:
	722	assert allclose(domain.quantities[name].explicit_update, 0)
	723	assert allclose(domain.quantities[name].semi_implicit_update, 0)
	724
	725	domain.compute_forcing_terms()
	726
	727	assert allclose(domain.quantities['level'].explicit_update, 0)
	728	assert allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
	729	assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
	730
	731	assert allclose(domain.quantities['level'].semi_implicit_update, 0)
	732	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
	733	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
	734
	735	#Create some momentum for friction to work with
	736	domain.set_quantity('xmomentum', 1)
	737	S = -g * eta2 / h(7.0/3)
	738
	739	domain.compute_forcing_terms()
	740	assert allclose(domain.quantities['level'].semi_implicit_update, 0)
	741	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
	742	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
	743
	744	#A more complex example
	745	domain.quantities['level'].semi_implicit_update[:] = 0.0
	746	domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
	747	domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
	748
	749	domain.set_quantity('xmomentum', 3)
	750	domain.set_quantity('ymomentum', 4)
	751
	752	S = -g * eta*2 5 / h**(7.0/3)
	753
	754
	755	domain.compute_forcing_terms()
	756
	757	assert allclose(domain.quantities['level'].semi_implicit_update, 0)
	758	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
	759	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)
	760
	761	def test_constant_wind_stress(self):
	762	from config import rho_a, rho_w, eta_w
	763	from math import pi, cos, sin, sqrt
	764
	765	a = [0.0, 0.0]
	766	b = [0.0, 2.0]
	767	c = [2.0, 0.0]
	768	d = [0.0, 4.0]
	769	e = [2.0, 2.0]
	770	f = [4.0, 0.0]
	771
	772	points = [a, b, c, d, e, f]
	773	#bac, bce, ecf, dbe
	774	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	775
	776
	777	domain = Domain(points, vertices)
	778
	779	#Flat surface with 1m of water
	780	domain.set_quantity('elevation', 0)
	781	domain.set_quantity('level', 1.0)
	782	domain.set_quantity('friction', 0)
	783
	784	Br = Reflective_boundary(domain)
	785	domain.set_boundary({'exterior': Br})
	786
	787	#Setup only one forcing term, constant wind stress
	788	s = 100
	789	phi = 135
	790	domain.forcing_terms = []
	791	domain.forcing_terms.append( Wind_stress(s, phi) )
	792
	793	domain.compute_forcing_terms()
	794
	795
	796	const = eta_w*rho_a/rho_w
	797
	798	#Convert to radians
	799	phi = phi*pi/180
	800
	801	#Compute velocity vector (u, v)
	802	u = s*cos(phi)
	803	v = s*sin(phi)
	804
	805	#Compute wind stress
	806	S = const * sqrt(u2 + v2)
	807
	808	assert allclose(domain.quantities['level'].explicit_update, 0)
	809	assert allclose(domain.quantities['xmomentum'].explicit_update, S*u)
	810	assert allclose(domain.quantities['ymomentum'].explicit_update, S*v)
	811
	812
	813	def test_variable_wind_stress(self):
	814	from config import rho_a, rho_w, eta_w
	815	from math import pi, cos, sin, sqrt
	816
	817	a = [0.0, 0.0]
	818	b = [0.0, 2.0]
	819	c = [2.0, 0.0]
	820	d = [0.0, 4.0]
	821	e = [2.0, 2.0]
	822	f = [4.0, 0.0]
	823
	824	points = [a, b, c, d, e, f]
	825	#bac, bce, ecf, dbe
	826	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	827
	828	domain = Domain(points, vertices)
	829
	830	#Flat surface with 1m of water
	831	domain.set_quantity('elevation', 0)
	832	domain.set_quantity('level', 1.0)
	833	domain.set_quantity('friction', 0)
	834
	835	Br = Reflective_boundary(domain)
	836	domain.set_boundary({'exterior': Br})
	837
	838
	839	domain.time = 5.54 #Take a random time (not zero)
	840
	841	#Setup only one forcing term, constant wind stress
	842	s = 100
	843	phi = 135
	844	domain.forcing_terms = []
	845	domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
	846
	847	domain.compute_forcing_terms()
	848
	849	#Compute reference solution
	850	const = eta_w*rho_a/rho_w
	851
	852	N = domain.number_of_elements
	853
	854	xc = domain.get_centroid_coordinates()
	855	t = domain.time
	856
[610]	857	x = xc[:,0]
	858	y = xc[:,1]
[623]	859	s_vec = speed(t,x,y)
	860	phi_vec = angle(t,x,y)
[610]	861
	862
[609]	863	for k in range(N):
	864	#Convert to radians
[610]	865	phi = phi_vec[k]*pi/180
	866	s = s_vec[k]
[609]	867
	868	#Compute velocity vector (u, v)
	869	u = s*cos(phi)
	870	v = s*sin(phi)
	871
	872	#Compute wind stress
	873	S = const * sqrt(u2 + v2)
	874
	875	assert allclose(domain.quantities['level'].explicit_update[k], 0)
	876	assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
	877	assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
	878
	879
	880
	881
[612]	882	def test_wind_stress_error_condition(self):
	883	"""Test that windstress reacts properly when forcing functions
	884	are wrong - e.g. returns a scalar
	885	"""
[609]	886
[612]	887	from config import rho_a, rho_w, eta_w
	888	from math import pi, cos, sin, sqrt
[609]	889
[612]	890	a = [0.0, 0.0]
	891	b = [0.0, 2.0]
	892	c = [2.0, 0.0]
	893	d = [0.0, 4.0]
	894	e = [2.0, 2.0]
	895	f = [4.0, 0.0]
	896
	897	points = [a, b, c, d, e, f]
	898	#bac, bce, ecf, dbe
	899	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	900
	901	domain = Domain(points, vertices)
	902
	903	#Flat surface with 1m of water
	904	domain.set_quantity('elevation', 0)
	905	domain.set_quantity('level', 1.0)
	906	domain.set_quantity('friction', 0)
	907
	908	Br = Reflective_boundary(domain)
	909	domain.set_boundary({'exterior': Br})
	910
	911
	912	domain.time = 5.54 #Take a random time (not zero)
	913
	914	#Setup only one forcing term, bad func
	915	domain.forcing_terms = []
	916
	917	try:
	918	domain.forcing_terms.append(Wind_stress(s = scalar_func,
	919	phi = angle))
	920	except AssertionError:
	921	pass
	922	else:
	923	msg = 'Should have raised exception'
	924	raise msg
[609]	925
	926
[612]	927	try:
	928	domain.forcing_terms.append(Wind_stress(s = speed,
	929	phi = scalar_func))
	930	except AssertionError:
	931	pass
	932	else:
	933	msg = 'Should have raised exception'
	934	raise msg
	935
	936	try:
	937	domain.forcing_terms.append(Wind_stress(s = speed,
	938	phi = 'xx'))
	939	except:
	940	pass
	941	else:
	942	msg = 'Should have raised exception'
	943	raise msg
	944
	945
[218]	946	#####################################################
	947	def test_first_order_extrapolator_const_z(self):
	948
	949	a = [0.0, 0.0]
	950	b = [0.0, 2.0]
	951	c = [2.0, 0.0]
	952	d = [0.0, 4.0]
	953	e = [2.0, 2.0]
	954	f = [4.0, 0.0]
	955
	956	points = [a, b, c, d, e, f]
	957	#bac, bce, ecf, dbe
	958	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	959
	960	domain = Domain(points, vertices)
	961	val0 = 2.+2.0/3
	962	val1 = 4.+4.0/3
	963	val2 = 8.+2.0/3
	964	val3 = 2.+8.0/3
	965
	966	zl=zr=-3.75 #Assume constant bed (must be less than level)
	967	domain.set_quantity('elevation', zl*ones( (4,3) ))
	968	domain.set_quantity('level', [[val0, val0-1, val0-2],
	969	[val1, val1+1, val1],
	970	[val2, val2-2, val2],
	971	[val3-0.5, val3, val3]])
	972
	973
	974
	975	domain.order = 1
	976	domain.distribute_to_vertices_and_edges()
	977
	978	#Check that centroid values were distributed to vertices
	979	C = domain.quantities['level'].centroid_values
	980	for i in range(3):
	981	assert allclose( domain.quantities['level'].vertex_values[:,i], C)
	982
	983
	984	def test_first_order_limiter_variable_z(self):
	985	#Check that first order limiter follows bed_slope
	986	from Numeric import alltrue, greater_equal
	987	from config import epsilon
	988
	989	a = [0.0, 0.0]
	990	b = [0.0, 2.0]
	991	c = [2.0,0.0]
	992	d = [0.0, 4.0]
	993	e = [2.0, 2.0]
	994	f = [4.0,0.0]
	995
	996	points = [a, b, c, d, e, f]
	997	#bac, bce, ecf, dbe
	998	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	999
	1000	domain = Domain(points, vertices)
	1001	val0 = 2.+2.0/3
	1002	val1 = 4.+4.0/3
	1003	val2 = 8.+2.0/3
	1004	val3 = 2.+8.0/3
	1005
	1006	domain.set_quantity('elevation', [[0,0,0], [6,0,0],
	1007	[6,6,6], [6,6,6]])
	1008	domain.set_quantity('level', [[val0, val0, val0],
	1009	[val1, val1, val1],
	1010	[val2, val2, val2],
	1011	[val3, val3, val3]])
	1012
	1013	E = domain.quantities['elevation'].vertex_values
	1014	L = domain.quantities['level'].vertex_values
	1015
	1016
	1017	#Check that some levels are not above elevation (within eps)
	1018	#- so that the limiter has something to work with
	1019	assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
	1020
	1021	domain.order = 1
	1022	domain.distribute_to_vertices_and_edges()
	1023
	1024	#Check that all levels are above elevation (within eps)
	1025	assert alltrue(alltrue(greater_equal(L,E-epsilon)))
	1026
	1027
	1028	#####################################################
	1029	def test_distribute_basic(self):
	1030	#Using test data generated by pyvolution-2
	1031	#Assuming no friction and flat bed (0.0)
	1032
	1033	a = [0.0, 0.0]
	1034	b = [0.0, 2.0]
	1035	c = [2.0, 0.0]
	1036	d = [0.0, 4.0]
	1037	e = [2.0, 2.0]
	1038	f = [4.0, 0.0]
	1039
	1040	points = [a, b, c, d, e, f]
	1041	#bac, bce, ecf, dbe
	1042	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1043
	1044	domain = Domain(points, vertices)
	1045
	1046	val0 = 2.
	1047	val1 = 4.
	1048	val2 = 8.
	1049	val3 = 2.
	1050
	1051	domain.set_quantity('level', [val0, val1, val2, val3], 'centroids')
	1052	L = domain.quantities['level'].vertex_values
	1053
	1054	#First order
	1055	domain.order = 1
	1056	domain.distribute_to_vertices_and_edges()
	1057	assert allclose(L[1], val1)
	1058
	1059	#Second order
	1060	domain.order = 2
	1061	domain.distribute_to_vertices_and_edges()
	1062	assert allclose(L[1], [2.2, 4.9, 4.9])
	1063
	1064
	1065
	1066	def test_distribute_away_from_bed(self):
	1067	#Using test data generated by pyvolution-2
	1068	#Assuming no friction and flat bed (0.0)
	1069
	1070	a = [0.0, 0.0]
	1071	b = [0.0, 2.0]
	1072	c = [2.0, 0.0]
	1073	d = [0.0, 4.0]
	1074	e = [2.0, 2.0]
	1075	f = [4.0, 0.0]
	1076
	1077	points = [a, b, c, d, e, f]
	1078	#bac, bce, ecf, dbe
	1079	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1080
	1081	domain = Domain(points, vertices)
	1082	L = domain.quantities['level'].vertex_values
	1083
	1084	def level(x,y):
	1085	return x**2
	1086
	1087	domain.set_quantity('level', level, 'centroids')
	1088
	1089	a, b = domain.quantities['level'].compute_gradients()
	1090	assert allclose(a[1], 3.33333334)
	1091	assert allclose(b[1], 0.0)
	1092
	1093	domain.order = 1
	1094	domain.distribute_to_vertices_and_edges()
	1095	assert allclose(L[1], 1.77777778)
	1096
	1097	domain.order = 2
	1098	domain.distribute_to_vertices_and_edges()
	1099	assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
	1100
	1101
	1102
	1103	def test_distribute_away_from_bed1(self):
	1104	#Using test data generated by pyvolution-2
	1105	#Assuming no friction and flat bed (0.0)
	1106
	1107	a = [0.0, 0.0]
	1108	b = [0.0, 2.0]
	1109	c = [2.0, 0.0]
	1110	d = [0.0, 4.0]
	1111	e = [2.0, 2.0]
	1112	f = [4.0, 0.0]
	1113
	1114	points = [a, b, c, d, e, f]
	1115	#bac, bce, ecf, dbe
	1116	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1117
	1118	domain = Domain(points, vertices)
	1119	L = domain.quantities['level'].vertex_values
	1120
	1121	def level(x,y):
	1122	return x4+y2
	1123
	1124	domain.set_quantity('level', level, 'centroids')
	1125	#print domain.quantities['level'].centroid_values
	1126
	1127	a, b = domain.quantities['level'].compute_gradients()
	1128	assert allclose(a[1], 25.18518519)
	1129	assert allclose(b[1], 3.33333333)
	1130
	1131	domain.order = 1
	1132	domain.distribute_to_vertices_and_edges()
	1133	assert allclose(L[1], 4.9382716)
	1134
	1135	domain.order = 2
	1136	domain.distribute_to_vertices_and_edges()
	1137	assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
	1138
	1139
	1140
	1141	def test_distribute_near_bed(self):
	1142	#Using test data generated by pyvolution-2
	1143	#Assuming no friction and flat bed (0.0)
	1144
	1145	a = [0.0, 0.0]
	1146	b = [0.0, 2.0]
	1147	c = [2.0, 0.0]
	1148	d = [0.0, 4.0]
	1149	e = [2.0, 2.0]
	1150	f = [4.0, 0.0]
	1151
	1152	points = [a, b, c, d, e, f]
	1153	#bac, bce, ecf, dbe
	1154	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1155
	1156	domain = Domain(points, vertices)
	1157
	1158
	1159	#Set up for a gradient of (3,0) at mid triangle
	1160	def slope(x, y):
	1161	return 10*x
	1162
	1163	h = 0.1
	1164	def level(x,y):
	1165	return slope(x,y)+h
	1166
	1167	domain.set_quantity('elevation', slope)
	1168	domain.set_quantity('level', level, 'centroids')
	1169
	1170	#print domain.quantities['elevation'].centroid_values
	1171	#print domain.quantities['level'].centroid_values
	1172
	1173	E = domain.quantities['elevation'].vertex_values
	1174	L = domain.quantities['level'].vertex_values
	1175
	1176	#print E
	1177	domain.order = 1
	1178	domain.distribute_to_vertices_and_edges()
[272]	1179	##assert allclose(L[1], [0.19999999, 20.05, 20.05])
	1180	assert allclose(L[1], [0.1, 20.1, 20.1])
	1181
[218]	1182	domain.order = 2
	1183	domain.distribute_to_vertices_and_edges()
	1184	assert allclose(L[1], [0.1, 20.1, 20.1])
	1185
	1186	def test_distribute_near_bed1(self):
	1187	#Using test data generated by pyvolution-2
	1188	#Assuming no friction and flat bed (0.0)
	1189
	1190	a = [0.0, 0.0]
	1191	b = [0.0, 2.0]
	1192	c = [2.0, 0.0]
	1193	d = [0.0, 4.0]
	1194	e = [2.0, 2.0]
	1195	f = [4.0, 0.0]
	1196
	1197	points = [a, b, c, d, e, f]
	1198	#bac, bce, ecf, dbe
	1199	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1200
	1201	domain = Domain(points, vertices)
	1202
	1203
	1204	#Set up for a gradient of (3,0) at mid triangle
	1205	def slope(x, y):
	1206	return x4+y2
	1207
	1208	h = 0.1
	1209	def level(x,y):
	1210	return slope(x,y)+h
	1211
	1212	domain.set_quantity('elevation', slope)
	1213	domain.set_quantity('level', level)
	1214
	1215	#print domain.quantities['elevation'].centroid_values
	1216	#print domain.quantities['level'].centroid_values
	1217
	1218	E = domain.quantities['elevation'].vertex_values
	1219	L = domain.quantities['level'].vertex_values
	1220
	1221	#print E
[272]	1222	domain.order = 1
[218]	1223	domain.distribute_to_vertices_and_edges()
[272]	1224	##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
	1225	assert allclose(L[1], [4.1, 16.1, 20.1])
[218]	1226
	1227	domain.order = 2
	1228	domain.distribute_to_vertices_and_edges()
	1229	assert allclose(L[1], [4.1, 16.1, 20.1])
	1230
[229]	1231
	1232
	1233	def test_second_order_distribute_real_data(self):
	1234	#Using test data generated by pyvolution-2
	1235	#Assuming no friction and flat bed (0.0)
	1236
	1237	a = [0.0, 0.0]
	1238	b = [0.0, 1.0/5]
	1239	c = [0.0, 2.0/5]
	1240	d = [1.0/5, 0.0]
	1241	e = [1.0/5, 1.0/5]
	1242	f = [1.0/5, 2.0/5]
	1243	g = [2.0/5, 2.0/5]
	1244
	1245	points = [a, b, c, d, e, f, g]
	1246	#bae, efb, cbf, feg
	1247	vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
[218]	1248
[229]	1249	domain = Domain(points, vertices)
	1250
	1251	def slope(x, y):
	1252	return -x/3
	1253
	1254	domain.set_quantity('elevation', slope)
	1255	domain.set_quantity('level',
	1256	[0.01298164, 0.00365611, 0.01440365, -0.0381856437096],
	1257	'centroids')
	1258	domain.set_quantity('xmomentum',
	1259	[0.00670439, 0.01263789, 0.00647805, 0.0178180740668],
	1260	'centroids')
	1261	domain.set_quantity('ymomentum',
	1262	[-7.23510980e-004, -6.30413883e-005, 6.30413883e-005, 0.000200907255866],
	1263	'centroids')
	1264
	1265	E = domain.quantities['elevation'].vertex_values
	1266	L = domain.quantities['level'].vertex_values
	1267	X = domain.quantities['xmomentum'].vertex_values
	1268	Y = domain.quantities['ymomentum'].vertex_values
	1269
	1270	#print E
	1271	domain.order = 2
	1272	domain.distribute_to_vertices_and_edges()
	1273
	1274	#print L[1,:]
	1275	#print X[1,:]
	1276	#print Y[1,:]
	1277
	1278	assert allclose(L[1,:], [-0.00825735775384, -0.00801881482869, 0.0272445025825])
	1279	assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717])
	1280	assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
	1281
	1282
[218]	1283	def test_second_order_flat_bed_onestep(self):
	1284
	1285	from mesh_factory import rectangular
	1286	from shallow_water import Domain, Reflective_boundary,\
	1287	Dirichlet_boundary, Constant_height
	1288	from Numeric import array
	1289
	1290	#Create basic mesh
	1291	points, vertices, boundary = rectangular(6, 6)
	1292
	1293	#Create shallow water domain
	1294	domain = Domain(points, vertices, boundary)
	1295	domain.smooth = False
	1296	domain.default_order=2
	1297
	1298	# Boundary conditions
	1299	Br = Reflective_boundary(domain)
	1300	Bd = Dirichlet_boundary([0.1, 0., 0.])
	1301	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
	1302
	1303	domain.check_integrity()
	1304
	1305	#Evolution
	1306	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
	1307	pass# domain.write_time()
	1308
	1309	#Data from earlier version of pyvolution
	1310	assert allclose(domain.min_timestep, 0.0396825396825)
	1311	assert allclose(domain.max_timestep, 0.0396825396825)
	1312
	1313	assert allclose(domain.quantities['level'].centroid_values[:12],
	1314	[0.00171396, 0.02656103, 0.00241523, 0.02656103,
	1315	0.00241523, 0.02656103, 0.00241523, 0.02656103,
	1316	0.00241523, 0.02656103, 0.00241523, 0.0272623])
	1317
	1318	domain.distribute_to_vertices_and_edges()
	1319	assert allclose(domain.quantities['level'].vertex_values[:12,0],
	1320	[0.0001714, 0.02656103, 0.00024152,
	1321	0.02656103, 0.00024152, 0.02656103,
	1322	0.00024152, 0.02656103, 0.00024152,
	1323	0.02656103, 0.00024152, 0.0272623])
	1324
	1325	assert allclose(domain.quantities['level'].vertex_values[:12,1],
	1326	[0.00315012, 0.02656103, 0.00024152, 0.02656103,
	1327	0.00024152, 0.02656103, 0.00024152, 0.02656103,
	1328	0.00024152, 0.02656103, 0.00040506, 0.0272623])
	1329
	1330	assert allclose(domain.quantities['level'].vertex_values[:12,2],
	1331	[0.00182037, 0.02656103, 0.00676264,
	1332	0.02656103, 0.00676264, 0.02656103,
	1333	0.00676264, 0.02656103, 0.00676264,
	1334	0.02656103, 0.0065991, 0.0272623])
	1335
	1336	assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
	1337	[0.00113961, 0.01302432, 0.00148672,
	1338	0.01302432, 0.00148672, 0.01302432,
	1339	0.00148672, 0.01302432, 0.00148672 ,
	1340	0.01302432, 0.00148672, 0.01337143])
	1341
	1342	assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
	1343	[-2.91240050e-004 , 1.22721531e-004,
	1344	-1.22721531e-004, 1.22721531e-004 ,
	1345	-1.22721531e-004, 1.22721531e-004,
	1346	-1.22721531e-004 , 1.22721531e-004,
	1347	-1.22721531e-004, 1.22721531e-004,
	1348	-1.22721531e-004, -4.57969873e-005])
	1349
	1350
	1351
	1352	def test_second_order_flat_bed_moresteps(self):
	1353
	1354	from mesh_factory import rectangular
	1355	from shallow_water import Domain, Reflective_boundary,\
	1356	Dirichlet_boundary, Constant_height
	1357	from Numeric import array
	1358
	1359	#Create basic mesh
	1360	points, vertices, boundary = rectangular(6, 6)
	1361
	1362	#Create shallow water domain
	1363	domain = Domain(points, vertices, boundary)
	1364	domain.smooth = False
	1365	domain.default_order=2
	1366
	1367	# Boundary conditions
	1368	Br = Reflective_boundary(domain)
	1369	Bd = Dirichlet_boundary([0.1, 0., 0.])
	1370	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
	1371
	1372	domain.check_integrity()
	1373
	1374	#Evolution
	1375	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
	1376	pass
	1377
	1378	#Data from earlier version of pyvolution
	1379	#assert allclose(domain.min_timestep, 0.0396825396825)
	1380	#assert allclose(domain.max_timestep, 0.0396825396825)
	1381	#print domain.quantities['level'].centroid_values
	1382
[221]	1383
	1384	def test_flatbed_first_order(self):
	1385	from mesh_factory import rectangular
	1386	from shallow_water import Domain,\
	1387	Reflective_boundary, Dirichlet_boundary
	1388
	1389	from Numeric import array
	1390
	1391	#Create basic mesh
	1392	N = 8
	1393	points, vertices, boundary = rectangular(N, N)
	1394
	1395	#Create shallow water domain
	1396	domain = Domain(points, vertices, boundary)
	1397	domain.smooth = False
	1398	domain.visualise = False
	1399	domain.default_order=1
	1400
	1401	# Boundary conditions
	1402	Br = Reflective_boundary(domain)
	1403	Bd = Dirichlet_boundary([0.2,0.,0.])
[218]	1404
[221]	1405	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
	1406	domain.check_integrity()
	1407
	1408
	1409	#Evolution
	1410	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
	1411	pass
	1412	#domain.write_time()
	1413
[443]	1414	#FIXME: These numbers were from version before 25/10
	1415	#assert allclose(domain.min_timestep, 0.0140413643926)
	1416	#assert allclose(domain.max_timestep, 0.0140947355753)
[221]	1417
	1418	for i in range(3):
[443]	1419	#assert allclose(domain.quantities['level'].edge_values[:4,i],
	1420	# [0.10730244,0.12337617,0.11200126,0.12605666])
[221]	1421
	1422	assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
	1423	[0.07610894,0.06901572,0.07284461,0.06819712])
	1424
[443]	1425	#assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
	1426	# [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
[218]	1427
[221]	1428	def test_flatbed_second_order(self):
	1429	from mesh_factory import rectangular
	1430	from shallow_water import Domain,\
	1431	Reflective_boundary, Dirichlet_boundary
	1432
	1433	from Numeric import array
	1434
	1435	#Create basic mesh
	1436	N = 8
	1437	points, vertices, boundary = rectangular(N, N)
	1438
	1439	#Create shallow water domain
	1440	domain = Domain(points, vertices, boundary)
	1441	domain.smooth = False
	1442	domain.visualise = False
	1443	domain.default_order=2
[272]	1444	#domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
[221]	1445
	1446	# Boundary conditions
	1447	Br = Reflective_boundary(domain)
	1448	Bd = Dirichlet_boundary([0.2,0.,0.])
	1449
	1450	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
	1451	domain.check_integrity()
	1452
	1453	#Evolution
	1454	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
	1455	pass
	1456
	1457
	1458	assert allclose(domain.min_timestep, 0.0210448446782)
	1459	assert allclose(domain.max_timestep, 0.0210448446782)
	1460
	1461	#print domain.quantities['level'].vertex_values[:4,0]
	1462	#print domain.quantities['xmomentum'].vertex_values[:4,0]
	1463	#print domain.quantities['ymomentum'].vertex_values[:4,0]
[443]	1464
	1465	#FIXME: These numbers were from version before 25/10
	1466	#assert allclose(domain.quantities['level'].vertex_values[:4,0],
	1467	# [0.00101913,0.05352143,0.00104852,0.05354394])
[221]	1468
	1469	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
[272]	1470	[ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
	1471
	1472	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
	1473	# [0.00090581,0.03685719,0.00088303,0.03687313])
	1474
[221]	1475	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
	1476	[-0.00139463,0.0006156,-0.00060364,0.00061827])
	1477
	1478
	1479
	1480	def test_flatbed_second_order_distribute(self):
	1481	#Use real data from pyvolution 2
	1482	#painfully setup and extracted.
	1483	from mesh_factory import rectangular
	1484	from shallow_water import Domain,\
	1485	Reflective_boundary, Dirichlet_boundary
	1486
	1487	from Numeric import array
	1488
	1489	#Create basic mesh
	1490	N = 8
	1491	points, vertices, boundary = rectangular(N, N)
	1492
	1493	#Create shallow water domain
	1494	domain = Domain(points, vertices, boundary)
	1495	domain.smooth = False
	1496	domain.visualise = False
	1497	domain.default_order=domain.order=2
	1498
	1499	# Boundary conditions
	1500	Br = Reflective_boundary(domain)
	1501	Bd = Dirichlet_boundary([0.2,0.,0.])
	1502
	1503	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
	1504	domain.check_integrity()
	1505
	1506
	1507
	1508	for V in [False, True]:
	1509	if V:
	1510	#Set centroids as if system had been evolved
	1511	L = zeros(2NN, Float)
	1512	L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
	1513	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
	1514	1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
	1515	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
	1516	1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
	1517	5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
	1518	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
	1519	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
	1520	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
	1521	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
	1522	0.00000000e+000, 5.57305948e-005]
	1523
	1524	X = zeros(2NN, Float)
	1525	X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
	1526	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
	1527	8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
	1528	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
	1529	8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
	1530	3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
	1531	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
	1532	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
	1533	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
	1534	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
	1535	0.00000000e+000, 4.57662812e-005]
	1536
	1537	Y = zeros(2NN, Float)
	1538	Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
	1539	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
	1540	-6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
	1541	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
	1542	-6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
	1543	-1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
	1544	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
	1545	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
	1546	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
	1547	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
	1548	0.00000000e+000, -2.57635178e-005]
	1549
	1550
	1551	domain.set_quantity('level', L, 'centroids')
	1552	domain.set_quantity('xmomentum', X, 'centroids')
	1553	domain.set_quantity('ymomentum', Y, 'centroids')
	1554
	1555	domain.check_integrity()
	1556	else:
	1557	#Evolution
	1558	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
	1559	pass
	1560	assert allclose(domain.min_timestep, 0.0210448446782)
	1561	assert allclose(domain.max_timestep, 0.0210448446782)
	1562
	1563
	1564	#Centroids were correct but not vertices.
	1565	#Hence the check of distribute below.
	1566	assert allclose(domain.quantities['level'].centroid_values[:4],
	1567	[0.00721206,0.05352143,0.01009108,0.05354394])
	1568
	1569	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
	1570	[0.00648352,0.03685719,0.00850733,0.03687313])
	1571
	1572	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
	1573	[-0.00139463,0.0006156,-0.00060364,0.00061827])
	1574
	1575	#print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
	1576	#print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
	1577
[272]	1578	#assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
	1579	##print domain.quantities['xmomentum'].centroid_values[17], V
	1580	##print
	1581	if not V:
	1582	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
	1583	else:
	1584	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
[221]	1585
	1586	import copy
	1587	XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
	1588	assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
	1589
	1590	domain.distribute_to_vertices_and_edges()
	1591
[272]	1592	#assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
[221]	1593
	1594	assert allclose(domain.quantities['xmomentum'].centroid_values[17],
[272]	1595	0.0)
[221]	1596
	1597
[443]	1598	#FIXME: These numbers were from version before 25/10
	1599	#assert allclose(domain.quantities['level'].vertex_values[:4,0],
	1600	# [0.00101913,0.05352143,0.00104852,0.05354394])
[221]	1601
	1602	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
	1603	[-0.00139463,0.0006156,-0.00060364,0.00061827])
	1604
[272]	1605
	1606	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
	1607	[0.00064835,0.03685719,0.00085073,0.03687313])
	1608
	1609
	1610	#NB NO longer relvant:
	1611
[221]	1612	#This was the culprit. First triangles vertex 0 had an
	1613	#x-momentum of 0.0064835 instead of 0.00090581 and
	1614	#third triangle had 0.00850733 instead of 0.00088303
	1615	#print domain.quantities['xmomentum'].vertex_values[:4,0]
	1616
[272]	1617	#print domain.quantities['xmomentum'].vertex_values[:4,0]
	1618	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
	1619	# [0.00090581,0.03685719,0.00088303,0.03687313])
[221]	1620
	1621
[272]	1622
[221]	1623
	1624
[218]	1625	def test_bedslope_problem_first_order(self):
	1626
	1627	from mesh_factory import rectangular
	1628	from shallow_water import Domain, Reflective_boundary, Constant_height
	1629	from Numeric import array
	1630
	1631	#Create basic mesh
	1632	points, vertices, boundary = rectangular(6, 6)
	1633
	1634	#Create shallow water domain
	1635	domain = Domain(points, vertices, boundary)
	1636	domain.smooth = False
	1637	domain.default_order=1
	1638
	1639	#Bed-slope and friction
	1640	def x_slope(x, y):
	1641	return -x/3
	1642
	1643	domain.set_quantity('elevation', x_slope)
	1644
	1645	# Boundary conditions
	1646	Br = Reflective_boundary(domain)
	1647	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1648
	1649	#Initial condition
	1650	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1651	domain.check_integrity()
	1652
	1653	#Evolution
	1654	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
	1655	pass# domain.write_time()
	1656
[272]	1657	assert allclose(domain.min_timestep, 0.050010003001)
	1658	assert allclose(domain.max_timestep, 0.050010003001)
[218]	1659
	1660
	1661	def test_bedslope_problem_first_order_moresteps(self):
	1662
	1663	from mesh_factory import rectangular
	1664	from shallow_water import Domain, Reflective_boundary, Constant_height
	1665	from Numeric import array
	1666
	1667	#Create basic mesh
	1668	points, vertices, boundary = rectangular(6, 6)
	1669
	1670	#Create shallow water domain
	1671	domain = Domain(points, vertices, boundary)
	1672	domain.smooth = False
	1673	domain.default_order=1
	1674
	1675	#Bed-slope and friction
	1676	def x_slope(x, y):
	1677	return -x/3
	1678
	1679	domain.set_quantity('elevation', x_slope)
	1680
	1681	# Boundary conditions
	1682	Br = Reflective_boundary(domain)
	1683	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1684
	1685	#Initial condition
	1686	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1687	domain.check_integrity()
	1688
	1689	#Evolution
	1690	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
	1691	pass# domain.write_time()
	1692
	1693	#Data from earlier version of pyvolution
[272]	1694	#print domain.quantities['level'].centroid_values
	1695
[218]	1696	assert allclose(domain.quantities['level'].centroid_values,
[272]	1697	[-0.02998628, -0.01520652, -0.03043492,
	1698	-0.0149132, -0.03004706, -0.01476251,
	1699	-0.0298215, -0.01467976, -0.02988158,
	1700	-0.01474662, -0.03036161, -0.01442995,
	1701	-0.07624583, -0.06297061, -0.07733792,
	1702	-0.06342237, -0.07695439, -0.06289595,
	1703	-0.07635559, -0.0626065, -0.07633628,
	1704	-0.06280072, -0.07739632, -0.06386738,
	1705	-0.12161738, -0.11028239, -0.1223796,
	1706	-0.11095953, -0.12189744, -0.11048616,
	1707	-0.12074535, -0.10987605, -0.12014311,
	1708	-0.10976691, -0.12096859, -0.11087692,
	1709	-0.16868259, -0.15868061, -0.16801135,
	1710	-0.1588003, -0.16674343, -0.15813323,
	1711	-0.16457595, -0.15693826, -0.16281096,
	1712	-0.15585154, -0.16283873, -0.15540068,
	1713	-0.17450362, -0.19919913, -0.18062882,
	1714	-0.19764131, -0.17783111, -0.19407213,
	1715	-0.1736915, -0.19053624, -0.17228678,
	1716	-0.19105634, -0.17920133, -0.1968828,
	1717	-0.14244395, -0.14604641, -0.14473537,
	1718	-0.1506107, -0.14510055, -0.14919522,
	1719	-0.14175896, -0.14560798, -0.13911658,
	1720	-0.14439383, -0.13924047, -0.14829043])
	1721
[218]	1722
[229]	1723	def test_bedslope_problem_second_order_one_step(self):
[218]	1724
	1725	from mesh_factory import rectangular
	1726	from shallow_water import Domain, Reflective_boundary, Constant_height
	1727	from Numeric import array
	1728
	1729	#Create basic mesh
	1730	points, vertices, boundary = rectangular(6, 6)
	1731
	1732	#Create shallow water domain
	1733	domain = Domain(points, vertices, boundary)
	1734	domain.smooth = False
	1735	domain.default_order=2
	1736
	1737	#Bed-slope and friction at vertices (and interpolated elsewhere)
	1738	def x_slope(x, y):
	1739	return -x/3
	1740
	1741	domain.set_quantity('elevation', x_slope)
	1742
	1743	# Boundary conditions
	1744	Br = Reflective_boundary(domain)
	1745	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1746
	1747	#Initial condition
	1748	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1749	domain.check_integrity()
	1750
	1751	assert allclose(domain.quantities['level'].centroid_values,
	1752	[0.01296296, 0.03148148, 0.01296296,
	1753	0.03148148, 0.01296296, 0.03148148,
	1754	0.01296296, 0.03148148, 0.01296296,
	1755	0.03148148, 0.01296296, 0.03148148,
	1756	-0.04259259, -0.02407407, -0.04259259,
	1757	-0.02407407, -0.04259259, -0.02407407,
	1758	-0.04259259, -0.02407407, -0.04259259,
	1759	-0.02407407, -0.04259259, -0.02407407,
	1760	-0.09814815, -0.07962963, -0.09814815,
	1761	-0.07962963, -0.09814815, -0.07962963,
	1762	-0.09814815, -0.07962963, -0.09814815,
	1763	-0.07962963, -0.09814815, -0.07962963,
	1764	-0.1537037 , -0.13518519, -0.1537037,
	1765	-0.13518519, -0.1537037, -0.13518519,
	1766	-0.1537037 , -0.13518519, -0.1537037,
	1767	-0.13518519, -0.1537037, -0.13518519,
	1768	-0.20925926, -0.19074074, -0.20925926,
	1769	-0.19074074, -0.20925926, -0.19074074,
	1770	-0.20925926, -0.19074074, -0.20925926,
	1771	-0.19074074, -0.20925926, -0.19074074,
	1772	-0.26481481, -0.2462963, -0.26481481,
	1773	-0.2462963, -0.26481481, -0.2462963,
	1774	-0.26481481, -0.2462963, -0.26481481,
	1775	-0.2462963, -0.26481481, -0.2462963])
	1776
	1777
	1778	#print domain.quantities['level'].extrapolate_second_order()
	1779	#domain.distribute_to_vertices_and_edges()
	1780	#print domain.quantities['level'].vertex_values[:,0]
	1781
	1782	#Evolution
[229]	1783	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
[218]	1784	#domain.write_time()
	1785	pass
	1786
	1787
[229]	1788	#print domain.quantities['level'].centroid_values
	1789	assert allclose(domain.quantities['level'].centroid_values,
	1790	[0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
	1791	0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
	1792	0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
	1793	-0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
	1794	-0.04186556, -0.0248011, -0.04186556, -0.02255593,
	1795	-0.09966629, -0.08035666, -0.09742112, -0.08035666,
	1796	-0.09742112, -0.08035666, -0.09742112, -0.08035666,
	1797	-0.09742112, -0.08035666, -0.09742112, -0.07811149,
	1798	-0.15522185, -0.13591222, -0.15297667, -0.13591222,
	1799	-0.15297667, -0.13591222, -0.15297667, -0.13591222,
	1800	-0.15297667, -0.13591222, -0.15297667, -0.13366704,
	1801	-0.2107774, -0.19146777, -0.20853223, -0.19146777,
	1802	-0.20853223, -0.19146777, -0.20853223, -0.19146777,
	1803	-0.20853223, -0.19146777, -0.20853223, -0.1892226,
	1804	-0.26120669, -0.24776246, -0.25865535, -0.24776246,
	1805	-0.25865535, -0.24776246, -0.25865535, -0.24776246,
	1806	-0.25865535, -0.24776246, -0.25865535, -0.24521113])
	1807
	1808
	1809
	1810	def test_bedslope_problem_second_order_two_steps(self):
	1811
	1812	from mesh_factory import rectangular
	1813	from shallow_water import Domain, Reflective_boundary, Constant_height
	1814	from Numeric import array
	1815
	1816	#Create basic mesh
	1817	points, vertices, boundary = rectangular(6, 6)
	1818
	1819	#Create shallow water domain
	1820	domain = Domain(points, vertices, boundary)
	1821	domain.smooth = False
	1822	domain.default_order=2
	1823
	1824	#Bed-slope and friction at vertices (and interpolated elsewhere)
	1825	def x_slope(x, y):
	1826	return -x/3
	1827
	1828	domain.set_quantity('elevation', x_slope)
	1829
	1830	# Boundary conditions
	1831	Br = Reflective_boundary(domain)
	1832	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1833
	1834	#Initial condition
	1835	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1836	domain.check_integrity()
	1837
	1838	assert allclose(domain.quantities['level'].centroid_values,
	1839	[0.01296296, 0.03148148, 0.01296296,
	1840	0.03148148, 0.01296296, 0.03148148,
	1841	0.01296296, 0.03148148, 0.01296296,
	1842	0.03148148, 0.01296296, 0.03148148,
	1843	-0.04259259, -0.02407407, -0.04259259,
	1844	-0.02407407, -0.04259259, -0.02407407,
	1845	-0.04259259, -0.02407407, -0.04259259,
	1846	-0.02407407, -0.04259259, -0.02407407,
	1847	-0.09814815, -0.07962963, -0.09814815,
	1848	-0.07962963, -0.09814815, -0.07962963,
	1849	-0.09814815, -0.07962963, -0.09814815,
	1850	-0.07962963, -0.09814815, -0.07962963,
	1851	-0.1537037 , -0.13518519, -0.1537037,
	1852	-0.13518519, -0.1537037, -0.13518519,
	1853	-0.1537037 , -0.13518519, -0.1537037,
	1854	-0.13518519, -0.1537037, -0.13518519,
	1855	-0.20925926, -0.19074074, -0.20925926,
	1856	-0.19074074, -0.20925926, -0.19074074,
	1857	-0.20925926, -0.19074074, -0.20925926,
	1858	-0.19074074, -0.20925926, -0.19074074,
	1859	-0.26481481, -0.2462963, -0.26481481,
	1860	-0.2462963, -0.26481481, -0.2462963,
	1861	-0.26481481, -0.2462963, -0.26481481,
	1862	-0.2462963, -0.26481481, -0.2462963])
	1863
	1864
	1865	#print domain.quantities['level'].extrapolate_second_order()
	1866	#domain.distribute_to_vertices_and_edges()
	1867	#print domain.quantities['level'].vertex_values[:,0]
[218]	1868
[229]	1869	#Evolution
	1870	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
	1871	pass
	1872
	1873
	1874	#Data from earlier version of pyvolution ft=0.1
[218]	1875	assert allclose(domain.min_timestep, 0.0376895634803)
	1876	assert allclose(domain.max_timestep, 0.0415635655309)
	1877
	1878
[229]	1879	assert allclose(domain.quantities['level'].centroid_values,
	1880	[0.00855788, 0.01575204, 0.00994606, 0.01717072,
	1881	0.01005985, 0.01716362, 0.01005985, 0.01716299,
	1882	0.01007098, 0.01736248, 0.01216452, 0.02026776,
	1883	-0.04462374, -0.02479045, -0.04199789, -0.0229465,
	1884	-0.04184033, -0.02295693, -0.04184013, -0.02295675,
	1885	-0.04184486, -0.0228168, -0.04028876, -0.02036486,
	1886	-0.10029444, -0.08170809, -0.09772846, -0.08021704,
	1887	-0.09760006, -0.08022143, -0.09759984, -0.08022124,
	1888	-0.09760261, -0.08008893, -0.09603914, -0.07758209,
	1889	-0.15584152, -0.13723138, -0.15327266, -0.13572906,
	1890	-0.15314427, -0.13573349, -0.15314405, -0.13573331,
	1891	-0.15314679, -0.13560104, -0.15158523, -0.13310701,
	1892	-0.21208605, -0.19283913, -0.20955631, -0.19134189,
	1893	-0.20942821, -0.19134598, -0.20942799, -0.1913458,
	1894	-0.20943005, -0.19120952, -0.20781177, -0.18869401,
	1895	-0.25384082, -0.2463294, -0.25047649, -0.24464654,
	1896	-0.25031159, -0.24464253, -0.25031112, -0.24464253,
	1897	-0.25031463, -0.24454764, -0.24885323, -0.24286438])
	1898
	1899
	1900
	1901
	1902	def test_bedslope_problem_second_order_two_yieldsteps(self):
	1903
	1904	from mesh_factory import rectangular
	1905	from shallow_water import Domain, Reflective_boundary, Constant_height
	1906	from Numeric import array
	1907
	1908	#Create basic mesh
	1909	points, vertices, boundary = rectangular(6, 6)
	1910
	1911	#Create shallow water domain
	1912	domain = Domain(points, vertices, boundary)
	1913	domain.smooth = False
	1914	domain.default_order=2
	1915
	1916	#Bed-slope and friction at vertices (and interpolated elsewhere)
	1917	def x_slope(x, y):
	1918	return -x/3
	1919
	1920	domain.set_quantity('elevation', x_slope)
	1921
	1922	# Boundary conditions
	1923	Br = Reflective_boundary(domain)
	1924	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1925
	1926	#Initial condition
	1927	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1928	domain.check_integrity()
	1929
	1930	assert allclose(domain.quantities['level'].centroid_values,
	1931	[0.01296296, 0.03148148, 0.01296296,
	1932	0.03148148, 0.01296296, 0.03148148,
	1933	0.01296296, 0.03148148, 0.01296296,
	1934	0.03148148, 0.01296296, 0.03148148,
	1935	-0.04259259, -0.02407407, -0.04259259,
	1936	-0.02407407, -0.04259259, -0.02407407,
	1937	-0.04259259, -0.02407407, -0.04259259,
	1938	-0.02407407, -0.04259259, -0.02407407,
	1939	-0.09814815, -0.07962963, -0.09814815,
	1940	-0.07962963, -0.09814815, -0.07962963,
	1941	-0.09814815, -0.07962963, -0.09814815,
	1942	-0.07962963, -0.09814815, -0.07962963,
	1943	-0.1537037 , -0.13518519, -0.1537037,
	1944	-0.13518519, -0.1537037, -0.13518519,
	1945	-0.1537037 , -0.13518519, -0.1537037,
	1946	-0.13518519, -0.1537037, -0.13518519,
	1947	-0.20925926, -0.19074074, -0.20925926,
	1948	-0.19074074, -0.20925926, -0.19074074,
	1949	-0.20925926, -0.19074074, -0.20925926,
	1950	-0.19074074, -0.20925926, -0.19074074,
	1951	-0.26481481, -0.2462963, -0.26481481,
	1952	-0.2462963, -0.26481481, -0.2462963,
	1953	-0.26481481, -0.2462963, -0.26481481,
	1954	-0.2462963, -0.26481481, -0.2462963])
	1955
	1956
	1957	#print domain.quantities['level'].extrapolate_second_order()
	1958	#domain.distribute_to_vertices_and_edges()
	1959	#print domain.quantities['level'].vertex_values[:,0]
[218]	1960
[229]	1961	#Evolution
	1962	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): #0.05??
	1963	#domain.write_time()
	1964	pass
	1965
	1966
	1967
	1968	assert allclose(domain.quantities['level'].centroid_values,
	1969	[0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
	1970	0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
	1971	0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
	1972	-0.0229465, -0.04184033, -0.02295693, -0.04184013,
	1973	-0.02295675, -0.04184486, -0.0228168, -0.04028876,
	1974	-0.02036486, -0.10029444, -0.08170809, -0.09772846,
	1975	-0.08021704, -0.09760006, -0.08022143, -0.09759984,
	1976	-0.08022124, -0.09760261, -0.08008893, -0.09603914,
	1977	-0.07758209, -0.15584152, -0.13723138, -0.15327266,
	1978	-0.13572906, -0.15314427, -0.13573349, -0.15314405,
	1979	-0.13573331, -0.15314679, -0.13560104, -0.15158523,
	1980	-0.13310701, -0.21208605, -0.19283913, -0.20955631,
	1981	-0.19134189, -0.20942821, -0.19134598, -0.20942799,
	1982	-0.1913458, -0.20943005, -0.19120952, -0.20781177,
	1983	-0.18869401, -0.25384082, -0.2463294, -0.25047649,
	1984	-0.24464654, -0.25031159, -0.24464253, -0.25031112,
	1985	-0.24464253, -0.25031463, -0.24454764, -0.24885323,
	1986	-0.24286438])
	1987
	1988
	1989
	1990	def test_bedslope_problem_second_order_more_steps(self):
	1991
	1992	from mesh_factory import rectangular
	1993	from shallow_water import Domain, Reflective_boundary, Constant_height
	1994	from Numeric import array
	1995
	1996	#Create basic mesh
	1997	points, vertices, boundary = rectangular(6, 6)
	1998
	1999	#Create shallow water domain
	2000	domain = Domain(points, vertices, boundary)
	2001	domain.smooth = False
	2002	domain.default_order=2
	2003
	2004	#Bed-slope and friction at vertices (and interpolated elsewhere)
	2005	def x_slope(x, y):
	2006	return -x/3
	2007
	2008	domain.set_quantity('elevation', x_slope)
	2009
	2010	# Boundary conditions
	2011	Br = Reflective_boundary(domain)
	2012	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	2013
	2014	#Initial condition
	2015	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	2016	domain.check_integrity()
	2017
	2018	assert allclose(domain.quantities['level'].centroid_values,
	2019	[0.01296296, 0.03148148, 0.01296296,
	2020	0.03148148, 0.01296296, 0.03148148,
	2021	0.01296296, 0.03148148, 0.01296296,
	2022	0.03148148, 0.01296296, 0.03148148,
	2023	-0.04259259, -0.02407407, -0.04259259,
	2024	-0.02407407, -0.04259259, -0.02407407,
	2025	-0.04259259, -0.02407407, -0.04259259,
	2026	-0.02407407, -0.04259259, -0.02407407,
	2027	-0.09814815, -0.07962963, -0.09814815,
	2028	-0.07962963, -0.09814815, -0.07962963,
	2029	-0.09814815, -0.07962963, -0.09814815,
	2030	-0.07962963, -0.09814815, -0.07962963,
	2031	-0.1537037 , -0.13518519, -0.1537037,
	2032	-0.13518519, -0.1537037, -0.13518519,
	2033	-0.1537037 , -0.13518519, -0.1537037,
	2034	-0.13518519, -0.1537037, -0.13518519,
	2035	-0.20925926, -0.19074074, -0.20925926,
	2036	-0.19074074, -0.20925926, -0.19074074,
	2037	-0.20925926, -0.19074074, -0.20925926,
	2038	-0.19074074, -0.20925926, -0.19074074,
	2039	-0.26481481, -0.2462963, -0.26481481,
	2040	-0.2462963, -0.26481481, -0.2462963,
	2041	-0.26481481, -0.2462963, -0.26481481,
	2042	-0.2462963, -0.26481481, -0.2462963])
	2043
	2044
	2045	#print domain.quantities['level'].extrapolate_second_order()
	2046	#domain.distribute_to_vertices_and_edges()
	2047	#print domain.quantities['level'].vertex_values[:,0]
	2048
	2049	#Evolution
	2050	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
	2051	pass
	2052
	2053
	2054	assert allclose(domain.quantities['level'].centroid_values,
	2055	[-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
	2056	-0.02901975, -0.0141361, -0.02898816, -0.01418135, -0.02961409, -0.01403487,
	2057	-0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
	2058	-0.07571145, -0.06235231, -0.0756817, -0.06245309, -0.07652292, -0.06289946,
	2059	-0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934, -0.1107174,
	2060	-0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
	2061	-0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493, -0.15697919,
	2062	-0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
	2063	-0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
	2064	-0.1833593, -0.19505871, -0.1818806, -0.19418042, -0.18586159, -0.19576946,
	2065	-0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
	2066	-0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
	2067
	2068	assert allclose(domain.quantities['xmomentum'].centroid_values,
	2069	[ 0.00831121, 0.00317948, 0.00731797, 0.00334939, 0.00764717, 0.00348053,
	2070	0.00788729, 0.00356522, 0.00780649, 0.00341919, 0.00693525, 0.00310375,
	2071	0.02166196, 0.01421475, 0.02017737, 0.01316839, 0.02037015, 0.01368659,
	2072	0.02106, 0.01399161, 0.02074514, 0.01354935, 0.01887407, 0.0123113,
	2073	0.03775083, 0.02855197, 0.03689337, 0.02759782, 0.03732848, 0.02812072,
	2074	0.03872545, 0.02913348, 0.03880939, 0.02803804, 0.03546499, 0.0260039,
	2075	0.0632131, 0.04730634, 0.0576324, 0.04592336, 0.05790921, 0.04690514,
	2076	0.05986467, 0.04871165, 0.06170068, 0.04811572, 0.05657041, 0.04416292,
	2077	0.08489642, 0.07188097, 0.07835261, 0.06843406, 0.07986412, 0.0698247,
	2078	0.08201071, 0.07216756, 0.08378418, 0.07273624, 0.080399, 0.06645841,
	2079	0.01631548, 0.04691608, 0.0206632, 0.044409, 0.02115518, 0.04560305,
	2080	0.02160608, 0.04663725, 0.02174734, 0.04795559, 0.02281427, 0.05667111])
	2081
	2082
	2083	assert allclose(domain.quantities['ymomentum'].centroid_values,
	2084	[ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
	2085	-9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
	2086	-1.99445058e-004, -3.28493467e-004, 6.68217349e-005, -8.42042805e-006,
	2087	5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
	2088	-8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
	2089	-2.70160173e-004, -5.35485454e-004, 2.60780997e-004, 3.12145471e-005,
	2090	5.16189608e-004, 1.07069062e-004, 9.29989252e-005, -3.71211119e-004,
	2091	1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
	2092	-4.74025708e-004, -6.94463009e-004, 6.15092843e-005, 2.22106820e-004,
	2093	-6.29589294e-004, 2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
	2094	-4.17914641e-004, 6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
	2095	-1.67424308e-003, -7.39485066e-004, -1.59966988e-003, 5.68262838e-005,
	2096	-1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
	2097	-1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
	2098	-2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
	2099	-1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
	2100	-6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
	2101	-4.50964850e-003, -3.06319963e-003, 6.08950810e-004, -4.79537921e-004])
	2102
	2103
	2104
	2105
	2106	def test_temp_play(self):
	2107
	2108	from mesh_factory import rectangular
	2109	from shallow_water import Domain, Reflective_boundary, Constant_height
	2110	from Numeric import array
	2111
	2112	#Create basic mesh
	2113	points, vertices, boundary = rectangular(5, 5)
	2114
	2115	#Create shallow water domain
	2116	domain = Domain(points, vertices, boundary)
	2117	domain.smooth = False
	2118	domain.default_order=2
	2119
	2120	#Bed-slope and friction at vertices (and interpolated elsewhere)
	2121	def x_slope(x, y):
	2122	return -x/3
	2123
	2124	domain.set_quantity('elevation', x_slope)
	2125
	2126	# Boundary conditions
	2127	Br = Reflective_boundary(domain)
	2128	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	2129
	2130	#Initial condition
	2131	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	2132	domain.check_integrity()
	2133
	2134	#Evolution
	2135	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
	2136	pass
	2137
	2138	assert allclose(domain.quantities['level'].centroid_values[:4],
	2139	[0.00206836, 0.01296714, 0.00363415, 0.01438924])
	2140	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
	2141	[0.01360154, 0.00671133, 0.01264578, 0.00648503])
	2142	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
	2143	[-1.19201077e-003, -7.23647546e-004,
	2144	-6.39083123e-005, 6.29815168e-005])
	2145
	2146
[232]	2147	def test_complex_bed(self):
	2148	#No friction is tested here
[229]	2149
[232]	2150	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
	2151	Transmissive_boundary, Time_boundary,\
	2152	Weir_simple as Weir, Constant_height
	2153
	2154	from mesh_factory import rectangular
	2155	from Numeric import array
[218]	2156
[232]	2157	N = 12
	2158	points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
	2159	origin=(-0.07, 0))
	2160
	2161
	2162	domain = Domain(points, vertices, boundary)
	2163	domain.smooth = False
	2164	domain.visualise = False
	2165	domain.default_order=2
	2166
	2167
	2168	inflow_stage = 0.1
	2169	Z = Weir(inflow_stage)
	2170	domain.set_quantity('elevation', Z)
	2171
	2172	Br = Reflective_boundary(domain)
	2173	Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
	2174	domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
	2175
	2176	domain.set_quantity('level', Constant_height(Z, 0.))
	2177
	2178	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
	2179	pass
	2180
[272]	2181
	2182	#print domain.quantities['level'].centroid_values
	2183
[443]	2184	#FIXME: These numbers were from version before 25/10
	2185	#assert allclose(domain.quantities['level'].centroid_values,
	2186	# [3.95822638e-002, 5.61022588e-002, 4.66437868e-002, 5.73081011e-002,
	2187	# 4.72394613e-002, 5.74684939e-002, 4.74309483e-002, 5.77458084e-002,
	2188	# 4.80628177e-002, 5.85656225e-002, 4.90498542e-002, 6.02609831e-002,
	2189	# 1.18470315e-002, 1.75136443e-002, 1.18035266e-002, 2.15565695e-002,
	2190	# 1.31620268e-002, 2.14351640e-002, 1.32351076e-002, 2.15450687e-002,
	2191	# 1.36414028e-002, 2.24274619e-002, 1.51689511e-002, 2.21789655e-002,
	2192	# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
	2193	# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
	2194	# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003, 2.73985295e-004,
	2195	# 1.86166683e-001, 8.74070369e-002, 1.86166712e-001, 8.74035875e-002,
	2196	# 6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
	2197	# 6.11666725e-002, 8.73846774e-002, 1.86166697e-001, 8.74171550e-002,
	2198	# -4.83333333e-002, 1.18333333e-001, -4.83333333e-002, 1.18333333e-001,
	2199	# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
	2200	# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002, 1.18333333e-001,
	2201	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
	2202	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
	2203	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
	2204	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
	2205	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
	2206	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
	2207	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
	2208	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
	2209	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
	2210	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
	2211	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
	2212	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
	2213	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
	2214	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
	2215	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
	2216	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
	2217	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
	2218	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
	2219	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
	2220	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
	2221	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
[232]	2222
[272]	2223
[232]	2224
	2225
[218]	2226	#-------------------------------------------------------------
	2227	if __name__ == "__main__":
	2228	suite = unittest.makeSuite(TestCase,'test')
[646]	2229	#suite = unittest.makeSuite(TestCase,'test_boundary_conditionsII')
[218]	2230	runner = unittest.TextTestRunner()
	2231	runner.run(suite)
	2232

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

Download in other formats: