Context Navigation

source: inundation/ga/storm_surge/pyvolution/test_quantity.py @ 546

Last change on this file since 546 was 546, checked in by duncan, 20 years ago
geting and setting values with a set of triangles - used with region tags
File size: 24.4 KB

Rev	Line
[261]	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
[476]	9	from Numeric import allclose, array, ones, Float
[261]	10
	11
	12	class TestCase(unittest.TestCase):
	13	def setUp(self):
	14	from domain import Domain
	15
	16	a = [0.0, 0.0]
	17	b = [0.0, 2.0]
	18	c = [2.0,0.0]
	19	d = [0.0, 4.0]
	20	e = [2.0, 2.0]
	21	f = [4.0,0.0]
	22
	23	points = [a, b, c, d, e, f]
	24
	25	#bac, bce, ecf, dbe
	26	elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
	27
[297]	28	self.mesh1 = Domain(points[:3], [elements[0]])
[261]	29	self.mesh1.check_integrity()
	30
	31	self.mesh4 = Domain(points, elements)
	32	self.mesh4.check_integrity()
	33
	34	def tearDown(self):
	35	pass
	36	#print " Tearing down"
	37
	38
	39	def test_creation(self):
	40
	41	quantity = Quantity(self.mesh1, [[1,2,3]])
	42	assert allclose(quantity.vertex_values, [[1.,2.,3.]])
	43
	44	try:
	45	quantity = Quantity()
	46	except:
	47	pass
	48	else:
	49	raise 'Should have raised empty quantity exception'
	50
	51
	52	try:
	53	quantity = Quantity([1,2,3])
	54	except AssertionError:
	55	pass
	56	except:
	57	raise 'Should have raised "mising mesh object" error'
	58
	59
	60	def test_creation_zeros(self):
	61
	62	quantity = Quantity(self.mesh1)
	63	assert allclose(quantity.vertex_values, [[0.,0.,0.]])
	64
	65
	66	quantity = Quantity(self.mesh4)
	67	assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
	68	[0.,0.,0.], [0.,0.,0.]])
	69
	70
	71	def test_interpolation(self):
	72	quantity = Quantity(self.mesh1, [[1,2,3]])
	73	assert allclose(quantity.centroid_values, [2.0]) #Centroid
	74
	75	assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
	76
	77
	78	def test_interpolation2(self):
[265]	79	quantity = Conserved_quantity(self.mesh4,
[261]	80	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	81	assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
	82
[265]	83
	84	quantity.extrapolate_second_order()
	85
	86	assert allclose(quantity.vertex_values, [[2., 2., 2.],
	87	[3.+2./3, 6.+2./3, 4.+2./3],
	88	[7.5, 0.5, 1.],
	89	[-5, -2.5, 7.5]])
	90
	91	#print quantity.edge_values
[261]	92	assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
	93	[5., 5., 5.],
	94	[4.5, 4.5, 0.],
	95	[3.0, -1.5, -1.5]])
	96
	97	def test_boundary_allocation(self):
	98	quantity = Conserved_quantity(self.mesh4,
	99	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	100
	101	assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
	102
	103
	104	def test_set_values(self):
	105	quantity = Quantity(self.mesh4)
	106
	107
	108	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
	109	location = 'vertices')
	110	assert allclose(quantity.vertex_values,
	111	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	112	assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
	113	assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
	114	[5., 5., 5.],
	115	[4.5, 4.5, 0.],
	116	[3.0, -1.5, -1.5]])
	117
	118
	119	#Test default
	120	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	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	#Test centroids
	130	quantity.set_values([1,2,3,4], location = 'centroids')
	131	assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
	132
	133	#Test edges
	134	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
	135	location = 'edges')
	136	assert allclose(quantity.edge_values,
	137	[[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
	138
	139	#Test exceptions
	140	try:
	141	quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
	142	location = 'bas kamel tuba')
	143	except:
	144	pass
	145
	146
	147	try:
	148	quantity.set_values([[1,2,3], [0,0,9]])
	149	except AssertionError:
	150	pass
	151	except:
	152	raise 'should have raised Assertionerror'
	153
	154
	155
	156	def test_set_values_const(self):
	157	quantity = Quantity(self.mesh4)
	158
	159	quantity.set_values(1.0, location = 'vertices')
	160	assert allclose(quantity.vertex_values,
	161	[[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
	162	assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
	163	assert allclose(quantity.edge_values, [[1, 1, 1],
	164	[1, 1, 1],
	165	[1, 1, 1],
	166	[1, 1, 1]])
	167
	168
	169	quantity.set_values(2.0, location = 'centroids')
	170	assert allclose(quantity.centroid_values, [2, 2, 2, 2])
	171
	172	quantity.set_values(3.0, location = 'edges')
	173	assert allclose(quantity.edge_values, [[3, 3, 3],
	174	[3, 3, 3],
	175	[3, 3, 3],
	176	[3, 3, 3]])
	177
	178
	179	def test_set_values_func(self):
	180	quantity = Quantity(self.mesh4)
	181
	182	def f(x, y):
	183	return x+y
	184
	185	quantity.set_values(f, location = 'vertices')
	186	assert allclose(quantity.vertex_values,
	187	[[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
	188	assert allclose(quantity.centroid_values,
	189	[4.0/3, 8.0/3, 10.0/3, 10.0/3])
	190	assert allclose(quantity.edge_values,
	191	[[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
	192
	193
	194	quantity.set_values(f, location = 'centroids')
	195	assert allclose(quantity.centroid_values,
	196	[4.0/3, 8.0/3, 10.0/3, 10.0/3])
	197
	198
[324]	199
	200
	201	def test_set_vertex_values(self):
	202	quantity = Quantity(self.mesh4)
	203
	204
	205	quantity.set_vertex_values([0,1,2,3,4,5])
	206
	207
	208	assert allclose(quantity.vertex_values,
	209	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
	210
	211	assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) #Centroid
	212
	213	assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
	214	[3., 2.5, 1.5],
	215	[3.5, 4.5, 3.],
	216	[2.5, 3.5, 2]])
	217
	218
	219
[344]	220
	221	def test_set_vertex_values_using_general_interface(self):
	222	quantity = Quantity(self.mesh4)
	223
	224
	225	quantity.set_values([0,1,2,3,4,5])
	226
	227
	228	assert allclose(quantity.vertex_values,
	229	[[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
	230
	231	assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) #Centroid
	232
	233	assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
	234	[3., 2.5, 1.5],
	235	[3.5, 4.5, 3.],
	236	[2.5, 3.5, 2]])
	237
	238
	239
	240
[261]	241	def test_gradient(self):
	242	quantity = Conserved_quantity(self.mesh4)
	243
	244	#Set up for a gradient of (3,0) at mid triangle
	245	quantity.set_values([2.0, 4.0, 8.0, 2.0],
	246	location = 'centroids')
[262]	247
	248
	249
[261]	250	#Gradients
	251	a, b = quantity.compute_gradients()
	252
	253	#gradient bewteen t0 and t1 is undefined as det==0
	254	assert a[0] == 0.0
	255	assert b[0] == 0.0
	256	#The others are OK
	257	for i in range(1,4):
	258	assert a[i] == 3.0
	259	assert b[i] == 0.0
	260
	261
	262	quantity.extrapolate_second_order()
[262]	263
	264	#print quantity.vertex_values
[261]	265	assert allclose(quantity.vertex_values, [[2., 2., 2.],
	266	[0., 6., 6.],
	267	[6., 6., 12.],
	268	[0., 0., 6.]])
	269
[262]	270
[261]	271	def test_second_order_extrapolation2(self):
	272	quantity = Conserved_quantity(self.mesh4)
	273
	274	#Set up for a gradient of (3,1), f(x) = 3x+y
	275	quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
	276	location = 'centroids')
	277
	278	#Gradients
	279	a, b = quantity.compute_gradients()
	280
	281	#gradient bewteen t0 and t1 is undefined as det==0
	282	assert a[0] == 0.0
	283	assert b[0] == 0.0
	284	#The others are OK
	285	for i in range(1,4):
	286	assert allclose(a[i], 3.0)
	287	assert allclose(b[i], 1.0)
	288
	289	quantity.extrapolate_second_order()
[265]	290
	291	#print quantity.vertex_values
[261]	292	assert allclose(quantity.vertex_values[1,0], 2.0)
	293	assert allclose(quantity.vertex_values[1,1], 6.0)
	294	assert allclose(quantity.vertex_values[1,2], 8.0)
	295
	296
	297
	298	def test_first_order_extrapolator(self):
	299	quantity = Conserved_quantity(self.mesh4)
	300
	301	#Test centroids
	302	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
	303	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
	304
	305	#Extrapolate
	306	quantity.extrapolate_first_order()
	307
	308	#Check vertices but not edge values
	309	assert allclose(quantity.vertex_values,
	310	[[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
	311
	312
	313	def test_second_order_extrapolator(self):
	314	quantity = Conserved_quantity(self.mesh4)
	315
	316	#Set up for a gradient of (3,0) at mid triangle
	317	quantity.set_values([2.0, 4.0, 8.0, 2.0],
	318	location = 'centroids')
	319
	320
	321
	322	quantity.extrapolate_second_order()
	323	quantity.limit()
	324
	325
	326	#Assert that central triangle is limited by neighbours
	327	assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
	328	assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
	329
	330	assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
	331	assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
	332
	333	assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
	334	assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
	335
	336
	337	#Assert that quantities are conserved
	338	from Numeric import sum
	339	for k in range(quantity.centroid_values.shape[0]):
	340	assert allclose (quantity.centroid_values[k],
	341	sum(quantity.vertex_values[k,:])/3)
	342
	343
	344
	345
	346
	347	def test_limiter(self):
	348	quantity = Conserved_quantity(self.mesh4)
	349
	350	#Create a deliberate overshoot (e.g. from gradient computation)
	351	quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
	352
	353
	354	#Limit
	355	quantity.limit()
	356
	357	#Assert that central triangle is limited by neighbours
	358	assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
	359	assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
	360
	361	assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
	362	assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
	363
	364	assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
	365	assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
	366
	367
	368
	369	#Assert that quantities are conserved
	370	from Numeric import sum
	371	for k in range(quantity.centroid_values.shape[0]):
	372	assert allclose (quantity.centroid_values[k],
	373	sum(quantity.vertex_values[k,:])/3)
	374
	375
	376
	377	def test_distribute_first_order(self):
	378	quantity = Conserved_quantity(self.mesh4)
	379
	380	#Test centroids
	381	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
	382	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
	383
	384
	385	#Extrapolate
	386	quantity.extrapolate_first_order()
	387
	388	#Interpolate
	389	quantity.interpolate_from_vertices_to_edges()
	390
	391	assert allclose(quantity.vertex_values,
	392	[[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
	393	assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
	394	[3,3,3], [4, 4, 4]])
	395
	396
	397
	398	def test_update_explicit(self):
	399	quantity = Conserved_quantity(self.mesh4)
	400
	401	#Test centroids
	402	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
	403	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
	404
	405	#Set explicit_update
	406	quantity.explicit_update = array( [1.,1.,1.,1.] )
	407
	408	#Update with given timestep
	409	quantity.update(0.1)
	410
	411	x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
	412	assert allclose( quantity.centroid_values, x)
	413
[476]	414	def test_update_semi_implicit(self):
	415	quantity = Conserved_quantity(self.mesh4)
[261]	416
[476]	417	#Test centroids
	418	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
	419	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
[261]	420
[476]	421	#Set semi implicit update
	422	quantity.semi_implicit_update = array([1.,1.,1.,1.])
[261]	423
[476]	424	#Update with given timestep
	425	timestep = 0.1
	426	quantity.update(timestep)
[261]	427
[476]	428	sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
	429	denom = ones(4, Float)-timestep*sem
	430
	431	x = array([1, 2, 3, 4])/denom
	432	assert allclose( quantity.centroid_values, x)
[261]	433
	434
[476]	435	def test_both_updates(self):
	436	quantity = Conserved_quantity(self.mesh4)
[261]	437
[476]	438	#Test centroids
	439	quantity.set_values([1.,2.,3.,4.], location = 'centroids')
	440	assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
[261]	441
[476]	442	#Set explicit_update
	443	quantity.explicit_update = array( [4.,3.,2.,1.] )
	444
	445	#Set semi implicit update
	446	quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
[261]	447
[476]	448	#Update with given timestep
	449	timestep = 0.1
	450	quantity.update(0.1)
[261]	451
[476]	452	sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
	453	denom = ones(4, Float)-timestep*sem
	454
	455	x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
	456	x /= denom
	457	assert allclose( quantity.centroid_values, x)
	458
[261]	459
[275]	460
	461
	462	#Test smoothing
	463	def test_smoothing(self):
	464
	465	from mesh_factory import rectangular
	466	from shallow_water import Domain, Transmissive_boundary
	467	from Numeric import zeros, Float
	468	from util import mean
	469
	470	#Create basic mesh
	471	points, vertices, boundary = rectangular(2, 2)
	472
	473	#Create shallow water domain
	474	domain = Domain(points, vertices, boundary)
	475	domain.default_order=2
	476	domain.reduction = mean
	477
[261]	478
[275]	479	#Set some field values
	480	domain.set_quantity('elevation', lambda x,y: x)
	481	domain.set_quantity('friction', 0.03)
	482
	483
	484	######################
	485	# Boundary conditions
	486	B = Transmissive_boundary(domain)
	487	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
	488
	489
	490	######################
	491	#Initial condition - with jumps
	492
	493	bed = domain.quantities['elevation'].vertex_values
	494	level = zeros(bed.shape, Float)
	495
	496	h = 0.03
	497	for i in range(level.shape[0]):
	498	if i % 2 == 0:
	499	level[i,:] = bed[i,:] + h
	500	else:
	501	level[i,:] = bed[i,:]
	502
	503	domain.set_quantity('level', level)
	504
	505	level = domain.quantities['level']
	506
	507	#Get smoothed level
	508	A, V = level.get_vertex_values(xy=False, smooth=True)
	509	Q = level.vertex_values
	510
[282]	511
[281]	512	assert A.shape[0] == 9
	513	assert V.shape[0] == 8
	514	assert V.shape[1] == 3
	515
[275]	516	#First four points
	517	assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
	518	assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
	519	assert allclose(A[2], Q[3,0])
	520	assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
	521
	522	#Center point
	523	assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
	524	Q[5,0] + Q[6,2] + Q[7,1])/6)
	525
	526
	527	#Check V
	528	assert allclose(V[0,:], [3,4,0])
	529	assert allclose(V[1,:], [1,0,4])
	530	assert allclose(V[2,:], [4,5,1])
	531	assert allclose(V[3,:], [2,1,5])
	532	assert allclose(V[4,:], [6,7,3])
	533	assert allclose(V[5,:], [4,3,7])
	534	assert allclose(V[6,:], [7,8,4])
	535	assert allclose(V[7,:], [5,4,8])
	536
	537	#Get smoothed level with XY
	538	X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=True)
	539
	540	assert allclose(A, A1)
	541	assert allclose(V, V1)
	542
	543	#Check XY
	544	assert allclose(X[4], 0.5)
	545	assert allclose(Y[4], 0.5)
	546
	547	assert allclose(X[7], 1.0)
	548	assert allclose(Y[7], 0.5)
	549
	550
	551
	552
	553	def test_vertex_values_no_smoothing(self):
	554
	555	from mesh_factory import rectangular
	556	from shallow_water import Domain, Transmissive_boundary
	557	from Numeric import zeros, Float
	558	from util import mean
	559
	560
	561	#Create basic mesh
	562	points, vertices, boundary = rectangular(2, 2)
	563
	564	#Create shallow water domain
	565	domain = Domain(points, vertices, boundary)
	566	domain.default_order=2
	567	domain.reduction = mean
	568
	569
	570	#Set some field values
	571	domain.set_quantity('elevation', lambda x,y: x)
	572	domain.set_quantity('friction', 0.03)
	573
[282]	574
	575	######################
	576	#Initial condition - with jumps
	577
	578	bed = domain.quantities['elevation'].vertex_values
	579	level = zeros(bed.shape, Float)
	580
	581	h = 0.03
	582	for i in range(level.shape[0]):
	583	if i % 2 == 0:
	584	level[i,:] = bed[i,:] + h
	585	else:
	586	level[i,:] = bed[i,:]
	587
	588	domain.set_quantity('level', level)
	589
[275]	590	#Get level
	591	level = domain.quantities['level']
	592	A, V = level.get_vertex_values(xy=False, smooth=False)
[282]	593	Q = level.vertex_values.flat
[275]	594
	595	for k in range(8):
	596	assert allclose(A[k], Q[k])
	597
[282]	598
	599	for k in range(8):
	600	assert V[k, 0] == 3*k
	601	assert V[k, 1] == 3*k+1
	602	assert V[k, 2] == 3*k+2
	603
	604
	605
[275]	606	X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=False)
	607
	608
	609	assert allclose(A, A1)
	610	assert allclose(V, V1)
	611
	612	#Check XY
[282]	613	assert allclose(X[1], 0.5)
	614	assert allclose(Y[1], 0.5)
	615	assert allclose(X[4], 0.0)
	616	assert allclose(Y[4], 0.0)
	617	assert allclose(X[12], 1.0)
	618	assert allclose(Y[12], 0.0)
[275]	619
	620
[459]	621
	622	def set_array_values_by_index(self):
[275]	623
[459]	624	from mesh_factory import rectangular
	625	from shallow_water import Domain
	626	from Numeric import zeros, Float
[275]	627
[459]	628	#Create basic mesh
	629	points, vertices, boundary = rectangular(1, 1)
	630
	631	#Create shallow water domain
	632	domain = Domain(points, vertices, boundary)
	633	#print "domain.number_of_elements ",domain.number_of_elements
	634	quantity = Quantity(domain,[[1,1,1],[2,2,2]])
	635	value = [7]
	636	indexes = [1]
	637	quantity.set_array_values_by_index(value,
	638	location = 'centroids',
	639	indexes = indexes)
	640	#print "quantity.centroid_values",quantity.centroid_values
	641
	642	assert allclose(quantity.centroid_values, [1,7])
	643
	644	quantity.set_array_values([15,20,25], indexes = indexes)
	645	assert allclose(quantity.centroid_values, [1,20])
	646
	647	quantity.set_array_values([15,20,25], indexes = indexes)
	648	assert allclose(quantity.centroid_values, [1,20])
	649
	650	def test_setting_some_vertex_values(self):
[517]	651	"""
	652	set values based on triangle lists.
	653	"""
[459]	654	from mesh_factory import rectangular
	655	from shallow_water import Domain
	656	from Numeric import zeros, Float
	657
	658	#Create basic mesh
	659	points, vertices, boundary = rectangular(1, 3)
	660
	661	#Create shallow water domain
	662	domain = Domain(points, vertices, boundary)
	663	#print "domain.number_of_elements ",domain.number_of_elements
	664	quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
	665	[4,4,4],[5,5,5],[6,6,6]])
	666	value = [7]
	667	indexes = [1]
[461]	668	quantity.set_values(value,
[459]	669	location = 'centroids',
	670	indexes = indexes)
	671	#print "quantity.centroid_values",quantity.centroid_values
	672	assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
	673
	674	value = [[15,20,25]]
[461]	675	quantity.set_values(value, indexes = indexes)
[459]	676	#print "1 quantity.vertex_values",quantity.vertex_values
	677	assert allclose(quantity.vertex_values[1], value[0])
	678
	679	values = [10,100,50]
[461]	680	quantity.set_values(values, indexes = [0,1,5])
[459]	681	#print "2 quantity.vertex_values",quantity.vertex_values
	682	assert allclose(quantity.vertex_values[0], [10,10,10])
	683	assert allclose(quantity.vertex_values[5], [50,50,50])
	684	quantity.interpolate()
	685	#print "quantity.centroid_values",quantity.centroid_values
	686	assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
	687
	688	values = [[31,30,29],[400,400,400],[1000,999,998]]
[461]	689	quantity.set_values(values, indexes = [3,3,5])
[459]	690	quantity.interpolate()
	691	assert allclose(quantity.centroid_values, [10,100,3,400,5,999])
[531]	692
	693
	694	def test_getting_some_vertex_values(self):
	695	"""
	696	get values based on triangle lists.
	697	"""
	698	from mesh_factory import rectangular
	699	from shallow_water import Domain
	700	from Numeric import zeros, Float
	701
	702	#Create basic mesh
	703	points, vertices, boundary = rectangular(1, 3)
	704
[546]	705	#print "points",points
	706	#print "vertices",vertices
	707	#print "boundary",boundary
	708
[531]	709	#Create shallow water domain
	710	domain = Domain(points, vertices, boundary)
	711	#print "domain.number_of_elements ",domain.number_of_elements
	712	quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
	713	[4,4,4],[5,5,5],[6,6,6]])
	714	value = [7]
	715	indexes = [1]
	716	quantity.set_values(value,
	717	location = 'centroids',
	718	indexes = indexes)
	719	#print "quantity.centroid_values",quantity.centroid_values
	720	#print "quantity.get_values(location = 'centroids') ",\
	721	# quantity.get_values(location = 'centroids')
	722	assert allclose(quantity.centroid_values,
	723	quantity.get_values(location = 'centroids'))
	724
	725
	726	value = [[15,20,25]]
	727	quantity.set_values(value, indexes = indexes)
	728	#print "1 quantity.vertex_values",quantity.vertex_values
	729	assert allclose(quantity.vertex_values, quantity.get_values())
	730
	731	assert allclose(quantity.edge_values,
	732	quantity.get_values(location = 'edges'))
	733
	734	# get a subset of elements
	735	subset = quantity.get_values(location='centroids', indexes=[0,5])
	736	answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
	737	assert allclose(subset, answer)
	738
	739
	740	subset = quantity.get_values(location='edges', indexes=[0,5])
	741	answer = [quantity.edge_values[0],quantity.edge_values[5]]
	742	#print "subset",subset
	743	#print "answer",answer
	744	assert allclose(subset, answer)
	745
	746	subset = quantity.get_values( indexes=[1,5])
	747	answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
	748	#print "subset",subset
	749	#print "answer",answer
	750	assert allclose(subset, answer)
	751
	752
[546]	753
[261]	754	#-------------------------------------------------------------
	755	if __name__ == "__main__":
	756	suite = unittest.makeSuite(TestCase,'test')
	757	runner = unittest.TextTestRunner()
	758	runner.run(suite)
	759
	760

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

Download in other formats: