Context Navigation

source: anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py @ 5571

Last change on this file since 5571 was 5288, checked in by ole, 17 years ago
Implemeted and tested ticket:175 - flow_through_cross_section. Also updated manual
File size: 54.1 KB

Line
1	#!/usr/bin/env python
2
3
4
5	#FIXME: Seperate the tests for mesh and general_mesh
6
7	#FIXME (Ole): Maxe this test independent of anything that inherits from General_mesh (namely shallow_water)
8
9	import unittest
10	from math import sqrt
11
12	from neighbour_mesh import *
13	from mesh_factory import rectangular
14	from anuga.config import epsilon
15	from Numeric import allclose, array, Int
16
17	from anuga.coordinate_transforms.geo_reference import Geo_reference
18	from anuga.utilities.polygon import is_inside_polygon
19	from anuga.utilities.numerical_tools import ensure_numeric
20
21	def distance(x, y):
22	return sqrt( sum( (array(x)-array(y))**2 ))
23
24	class Test_Mesh(unittest.TestCase):
25	def setUp(self):
26	pass
27
28	def tearDown(self):
29	pass
30
31	def test_triangle_inputs(self):
32	points = [[0.0, 0.0], [4.0, 0.0], [0.0, 3.0]]
33	vertices = [0,1,2] #Wrong
34
35	try:
36	mesh = Mesh(points, vertices)
37	except:
38	pass
39	else:
40	msg = 'Should have raised exception'
41	raise msg
42
43
44	def test_basic_triangle(self):
45
46	a = [0.0, 0.0]
47	b = [4.0, 0.0]
48	c = [0.0, 3.0]
49
50	points = [a, b, c]
51	vertices = [[0,1,2]]
52	mesh = Mesh(points, vertices)
53
54	#Centroid
55	centroid = mesh.centroid_coordinates[0]
56	assert centroid[0] == 4.0/3
57	assert centroid[1] == 1.0
58
59	#Area
60	assert mesh.areas[0] == 6.0,\
61	'Area was %f, should have been 6.0' %mesh.areas[0]
62
63	#Normals
64	normals = mesh.get_normals()
65	assert allclose(normals[0, 0:2], [3.0/5, 4.0/5])
66	assert allclose(normals[0, 2:4], [-1.0, 0.0])
67	assert allclose(normals[0, 4:6], [0.0, -1.0])
68
69	assert allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])
70	assert allclose(mesh.get_normal(0,1), [-1.0, 0.0])
71	assert allclose(mesh.get_normal(0,2), [0.0, -1.0])
72
73	#Edge lengths
74	assert allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])
75
76
77	#Vertex coordinates
78	#V = mesh.get_vertex_coordinates()
79	#assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0])
80
81
82	V = mesh.get_vertex_coordinates()
83	assert allclose(V, [ [0.0, 0.0],
84	[4.0, 0.0],
85	[0.0, 3.0] ])
86
87	V0 = mesh.get_vertex_coordinate(0, 0)
88	assert allclose(V0, [0.0, 0.0])
89
90	V1 = mesh.get_vertex_coordinate(0, 1)
91	assert allclose(V1, [4.0, 0.0])
92
93	V2 = mesh.get_vertex_coordinate(0, 2)
94	assert allclose(V2, [0.0, 3.0])
95
96
97	#General tests:
98
99	#Test that points are arranged in a counter clock wise order etc
100	mesh.check_integrity()
101
102
103	#Test that the centroid is located 2/3 of the way
104	#from each vertex to the midpoint of the opposite side
105
106	V = mesh.get_vertex_coordinates()
107	x0 = V[0, 0]; y0 = V[0, 1]
108	x1 = V[1, 0]; y1 = V[1, 1]
109	x2 = V[2, 0]; y2 = V[2, 1]
110	#x0 = V[0,0]
111	#y0 = V[0,1]
112	#x1 = V[0,2]
113	#y1 = V[0,3]
114	#x2 = V[0,4]
115	#y2 = V[0,5]
116
117	m0 = [(x1 + x2)/2, (y1 + y2)/2]
118	m1 = [(x0 + x2)/2, (y0 + y2)/2]
119	m2 = [(x1 + x0)/2, (y1 + y0)/2]
120
121	d0 = distance(centroid, [x0, y0])
122	d1 = distance(m0, [x0, y0])
123	assert d0 == 2*d1/3
124	#
125	d0 = distance(centroid, [x1, y1])
126	d1 = distance(m1, [x1, y1])
127	assert abs(d0 - 2d1/3) < epsilon, '%e, %e' %(d0, 2d1/3)
128
129	d0 = distance(centroid, [x2, y2])
130	d1 = distance(m2, [x2, y2])
131	assert abs(d0 - 2d1/3) < epsilon, '%e, %e' %(d0, 2d1/3)
132
133	#Radius
134	d0 = distance(centroid, m0)
135	assert d0 == 5.0/6
136
137	d1 = distance(centroid, m1)
138	assert d1 == sqrt(73.0/36)
139
140	d2 = distance(centroid, m2)
141	assert d2 == sqrt(13.0/9)
142
143	assert mesh.radii[0] == min(d0, d1, d2)
144	assert mesh.radii[0] == 5.0/6
145
146
147	#Let x be the centroid of triangle abc.
148	#Test that areas of the three triangles axc, cxb, and bxa are equal.
149	points = [a, b, c, centroid]
150	vertices = [[0,3,2], [2,3,1], [1,3,0]]
151	new_mesh = Mesh(points, vertices)
152
153	assert new_mesh.areas[0] == new_mesh.areas[1]
154	assert new_mesh.areas[1] == new_mesh.areas[2]
155	assert new_mesh.areas[1] == new_mesh.areas[2]
156
157	assert new_mesh.areas[1] == mesh.areas[0]/3
158
159
160
161	def test_general_triangle(self):
162	a = [2.0, 1.0]
163	b = [6.0, 2.0]
164	c = [1.0, 3.0]
165
166	points = [a, b, c]
167	vertices = [[0,1,2]]
168
169	mesh = Mesh(points, vertices)
170	centroid = mesh.centroid_coordinates[0]
171
172
173	#Test that the centroid is located 2/3 of the way
174	#from each vertex to the midpoint of the opposite side
175
176	V = mesh.get_vertex_coordinates()
177	x0 = V[0, 0]; y0 = V[0, 1]
178	x1 = V[1, 0]; y1 = V[1, 1]
179	x2 = V[2, 0]; y2 = V[2, 1]
180
181	#x0 = V[0,0]
182	#y0 = V[0,1]
183	#x1 = V[0,2]
184	#y1 = V[0,3]
185	#x2 = V[0,4]
186	#y2 = V[0,5]
187
188	m0 = [(x1 + x2)/2, (y1 + y2)/2]
189	m1 = [(x0 + x2)/2, (y0 + y2)/2]
190	m2 = [(x1 + x0)/2, (y1 + y0)/2]
191
192	d0 = distance(centroid, [x0, y0])
193	d1 = distance(m0, [x0, y0])
194	assert abs(d0 - 2d1/3) < epsilon, '%e, %e' %(d0, 2d1/3)
195	#
196	d0 = distance(centroid, [x1, y1])
197	d1 = distance(m1, [x1, y1])
198	assert abs(d0 - 2d1/3) < epsilon, '%e, %e' %(d0, 2d1/3)
199
200	d0 = distance(centroid, [x2, y2])
201	d1 = distance(m2, [x2, y2])
202	assert abs(d0 - 2d1/3) < epsilon, '%e, %e' %(d0, 2d1/3)
203
204	#Radius
205	d0 = distance(centroid, m0)
206	d1 = distance(centroid, m1)
207	d2 = distance(centroid, m2)
208	assert mesh.radii[0] == min(d0, d1, d2)
209
210
211
212	#Let x be the centroid of triangle abc.
213	#Test that areas of the three triangles axc, cxb, and bxa are equal.
214
215	points = [a, b, c, centroid]
216	vertices = [[0,3,2], [2,3,1], [1,3,0]]
217	new_mesh = Mesh(points, vertices)
218
219	assert new_mesh.areas[0] == new_mesh.areas[1]
220	assert new_mesh.areas[1] == new_mesh.areas[2]
221	assert new_mesh.areas[1] == new_mesh.areas[2]
222
223	assert new_mesh.areas[1] == mesh.areas[0]/3
224
225
226	#Test that points are arranged in a counter clock wise order
227	mesh.check_integrity()
228
229	def test_inscribed_circle_equilateral(self):
230	"""test that the radius is calculated correctly by mesh in the case of an equilateral triangle"""
231	a = [0.0, 0.0]
232	b = [2.0, 0.0]
233	c = [1.0, sqrt(3.0)]
234
235	points = [a, b, c]
236	vertices = [[0,1,2]]
237
238	mesh = Mesh(points, vertices,use_inscribed_circle=False)
239	assert allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work'
240
241	mesh = Mesh(points, vertices,use_inscribed_circle=True)
242	assert allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work'
243
244	def test_inscribed_circle_rightangle_triangle(self):
245	"""test that the radius is calculated correctly by mesh in the case of a right-angled triangle"""
246	a = [0.0, 0.0]
247	b = [4.0, 0.0]
248	c = [0.0, 3.0]
249
250	points = [a, b, c]
251	vertices = [[0,1,2]]
252
253	mesh = Mesh(points, vertices,use_inscribed_circle=False)
254	assert allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work'
255
256	mesh = Mesh(points, vertices,use_inscribed_circle=True)
257	assert allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work'
258
259
260	def test_two_triangles(self):
261	a = [0.0, 0.0]
262	b = [0.0, 2.0]
263	c = [2.0,0.0]
264	e = [2.0, 2.0]
265	points = [a, b, c, e]
266	vertices = [ [1,0,2], [1,2,3] ] #bac, bce
267	mesh = Mesh(points, vertices)
268
269	assert mesh.areas[0] == 2.0
270
271	assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
272
273
274	#Test that points are arranged in a counter clock wise order
275	mesh.check_integrity()
276
277
278
279	def test_more_triangles(self):
280
281	a = [0.0, 0.0]
282	b = [0.0, 2.0]
283	c = [2.0, 0.0]
284	d = [0.0, 4.0]
285	e = [2.0, 2.0]
286	f = [4.0, 0.0]
287
288	points = [a, b, c, d, e, f]
289	#bac, bce, ecf, dbe, daf, dae
290	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
291	mesh = Mesh(points, vertices)
292
293	#Test that points are arranged in a counter clock wise order
294	mesh.check_integrity()
295
296	assert mesh.areas[0] == 2.0
297	assert mesh.areas[1] == 2.0
298	assert mesh.areas[2] == 2.0
299	assert mesh.areas[3] == 2.0
300
301	assert mesh.edgelengths[1,0] == 2.0
302	assert mesh.edgelengths[1,1] == 2.0
303	assert mesh.edgelengths[1,2] == sqrt(8.0)
304
305	assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
306	assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3])
307	assert allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])
308	assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])
309
310	def test_mesh_and_neighbours(self):
311	a = [0.0, 0.0]
312	b = [0.0, 2.0]
313	c = [2.0,0.0]
314	d = [0.0, 4.0]
315	e = [2.0, 2.0]
316	f = [4.0,0.0]
317
318
319	points = [a, b, c, d, e, f]
320
321	#bac, bce, ecf, dbe
322	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
323	mesh = Mesh(points, vertices)
324
325	mesh.check_integrity()
326
327
328	T = mesh
329	tid = 0
330	assert T.number_of_boundaries[tid] == 2
331	assert T.neighbours[tid, 0] < 0 #Opposite point b (0,2)
332	assert T.neighbours[tid, 1] == 1 #Opposite point a (0,0)
333	assert T.neighbours[tid, 2] < 0 #Opposite point c (2,0)
334
335	tid = 1
336	assert T.number_of_boundaries[tid] == 0
337	assert T.neighbours[tid, 0] == 2 #Opposite point b (0,2)
338	assert T.neighbours[tid, 1] == 3 #Opposite point c (2,0)
339	assert T.neighbours[tid, 2] == 0 #Opposite point e (2,2)
340
341	tid = 2
342	assert T.number_of_boundaries[tid] == 2
343	assert T.neighbours[tid, 0] < 0 #Opposite point e (2,2)
344	assert T.neighbours[tid, 1] < 0 #Opposite point c (2,0)
345	assert T.neighbours[tid, 2] == 1 #Opposite point f (4,0)
346
347	tid = 3
348	assert T.number_of_boundaries[tid] == 2
349	assert T.neighbours[tid, 0] == 1 #Opposite point d (0,4)
350	assert T.neighbours[tid, 1] < 0 #Opposite point b (0,3)
351	assert T.neighbours[tid, 2] < 0 #Opposite point e (2,2)
352
353	#Neighbouring edges
354	tid = 0
355	assert T.neighbour_edges[tid, 0] < 0 #Opposite point b (0,2)
356	assert T.neighbour_edges[tid, 1] == 2 #Opposite point a (0,0)
357	assert T.neighbour_edges[tid, 2] < 0 #Opposite point c (2,0)
358
359	tid = 1
360	assert T.neighbour_edges[tid, 0] == 2 #Opposite point b (0,2)
361	assert T.neighbour_edges[tid, 1] == 0 #Opposite point c (2,0)
362	assert T.neighbour_edges[tid, 2] == 1 #Opposite point e (2,2)
363
364	tid = 2
365	assert T.neighbour_edges[tid, 0] < 0 #Opposite point e (2,2)
366	assert T.neighbour_edges[tid, 1] < 0 #Opposite point c (2,0)
367	assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0)
368
369	tid = 3
370	assert T.neighbour_edges[tid, 0] == 1 #Opposite point d (0,4)
371	assert T.neighbour_edges[tid, 1] < 0 #Opposite point b (0,3)
372	assert T.neighbour_edges[tid, 2] < 0 #Opposite point e (2,2)
373
374
375	def test_build_neighbour_structure_duplicates(self):
376	p0 = [-66.0, 14.0]
377	p1 = [14.0, -66.0]
378	p2 = [14.0, 14.0]
379	p3 = [60.0, 20.0]
380	p4 = [10.0, 60.0]
381	p5 = [60.0, 60.0]
382
383	points = [p0, p1, p2, p3, p4, p5]
384	triangles = [ [0, 1, 2],
385	[3, 2, 1],
386	[0, 2, 4],
387	[0, 2, 4],
388	[4, 2, 5],
389	[5, 2, 3]]
390	try:
391	mesh = Mesh(points, triangles)
392	except:
393	pass
394	else:
395	raise "triangle edge duplicates not caught"
396
397	def test_rectangular_mesh_basic(self):
398	M=1
399	N=1
400
401	points, vertices, boundary = rectangular(M, N)
402	mesh = Mesh(points, vertices, boundary)
403
404	#Test that points are arranged in a counter clock wise order
405	mesh.check_integrity()
406
407	M=2
408	N=2
409	points, vertices, boundary = rectangular(M, N)
410	mesh = Mesh(points, vertices, boundary)
411
412	#Test that points are arranged in a counter clock wise order
413	mesh.check_integrity()
414
415	#assert mesh.boundary[(7,1)] == 2 # top
416	assert mesh.boundary[(7,1)] == 'top' # top
417	assert mesh.boundary[(3,1)] == 'top' # top
418
419
420	def test_boundary_tags(self):
421
422
423	points, vertices, boundary = rectangular(4, 4)
424	mesh = Mesh(points, vertices, boundary)
425
426
427	#Test that points are arranged in a counter clock wise order
428	mesh.check_integrity()
429
430	#print mesh.get_boundary_tags()
431	#print mesh.boundary
432
433	for k in [1, 3, 5, 7]:
434	assert mesh.boundary[(k,2)] == 'left'
435
436	for k in [24, 26, 28, 30]:
437	assert mesh.boundary[(k,2)] == 'right'
438
439	for k in [7, 15, 23, 31]:
440	assert mesh.boundary[(k,1)] == 'top'
441	for k in [0, 8, 16, 24]:
442	assert mesh.boundary[(k,1)] == 'bottom'
443
444
445
446	def test_rectangular_mesh(self):
447	M=4
448	N=16
449	len1 = 100.0
450	len2 = 17.0
451
452	points, vertices, boundary = rectangular(M, N, len1, len2)
453	mesh = Mesh(points, vertices, boundary)
454
455	assert len(mesh) == 2MN
456
457	for i in range(len(mesh)):
458	assert mesh.areas[i] == len1len2/(2M*N)
459
460	hypo = sqrt((len1/M)2 + (len2/N)2) #hypothenuse
461	assert mesh.edgelengths[i, 0] == hypo
462	assert mesh.edgelengths[i, 1] == len1/M #x direction
463	assert mesh.edgelengths[i, 2] == len2/N #y direction
464
465	#Test that points are arranged in a counter clock wise order
466	mesh.check_integrity()
467
468
469	def test_rectangular_mesh2(self):
470	#Check that integers don't cause trouble
471	N = 16
472
473	points, vertices, boundary = rectangular(2*N, N, len1=10, len2=10)
474	mesh = Mesh(points, vertices, boundary)
475
476
477
478	def test_surrogate_neighbours(self):
479	a = [0.0, 0.0]
480	b = [0.0, 2.0]
481	c = [2.0,0.0]
482	d = [0.0, 4.0]
483	e = [2.0, 2.0]
484	f = [4.0,0.0]
485
486	points = [a, b, c, d, e, f]
487
488	#bac, bce, ecf, dbe
489	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
490	mesh = Mesh(points, vertices)
491	mesh.check_integrity()
492
493
494	T = mesh
495	tid = 0
496	assert T.number_of_boundaries[tid] == 2
497	assert T.surrogate_neighbours[tid, 0] == tid
498	assert T.surrogate_neighbours[tid, 1] == 1
499	assert T.surrogate_neighbours[tid, 2] == tid
500
501	tid = 1
502	assert T.number_of_boundaries[tid] == 0
503	assert T.surrogate_neighbours[tid, 0] == 2
504	assert T.surrogate_neighbours[tid, 1] == 3
505	assert T.surrogate_neighbours[tid, 2] == 0
506
507	tid = 2
508	assert T.number_of_boundaries[tid] == 2
509	assert T.surrogate_neighbours[tid, 0] == tid
510	assert T.surrogate_neighbours[tid, 1] == tid
511	assert T.surrogate_neighbours[tid, 2] == 1
512
513	tid = 3
514	assert T.number_of_boundaries[tid] == 2
515	assert T.surrogate_neighbours[tid, 0] == 1
516	assert T.surrogate_neighbours[tid, 1] == tid
517	assert T.surrogate_neighbours[tid, 2] == tid
518
519
520	def test_boundary_inputs(self):
521	a = [0.0, 0.0]
522	b = [0.0, 2.0]
523	c = [2.0,0.0]
524	d = [0.0, 4.0]
525	e = [2.0, 2.0]
526	f = [4.0,0.0]
527
528	points = [a, b, c, d, e, f]
529
530	#bac, bce, ecf, dbe
531	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
532
533	boundary = { (0, 0): 'First',
534	(0, 2): 'Second',
535	(2, 0): 'Third',
536	(2, 1): 'Fourth',
537	(3, 1): 'Fifth',
538	(3, 2): 'Sixth'}
539
540
541	mesh = Mesh(points, vertices, boundary)
542	mesh.check_integrity()
543
544
545	#Check enumeration
546	#for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
547	# b = -k-1
548	# assert mesh.neighbours[vol_id, edge_id] == b
549
550
551
552	def test_boundary_inputs_using_one_default(self):
553	a = [0.0, 0.0]
554	b = [0.0, 2.0]
555	c = [2.0,0.0]
556	d = [0.0, 4.0]
557	e = [2.0, 2.0]
558	f = [4.0,0.0]
559
560	points = [a, b, c, d, e, f]
561
562	#bac, bce, ecf, dbe
563	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
564
565	boundary = { (0, 0): 'First',
566	(0, 2): 'Second',
567	(2, 0): 'Third',
568	(2, 1): 'Fourth',
569	#(3, 1): 'Fifth', #Skip this
570	(3, 2): 'Sixth'}
571
572
573	mesh = Mesh(points, vertices, boundary)
574	mesh.check_integrity()
575
576	from anuga.config import default_boundary_tag
577	assert mesh.boundary[ (3, 1) ] == default_boundary_tag
578
579
580	#Check enumeration
581	#for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
582	# b = -k-1
583	# assert mesh.neighbours[vol_id, edge_id] == b
584
585	def test_boundary_inputs_using_all_defaults(self):
586	a = [0.0, 0.0]
587	b = [0.0, 2.0]
588	c = [2.0,0.0]
589	d = [0.0, 4.0]
590	e = [2.0, 2.0]
591	f = [4.0,0.0]
592
593	points = [a, b, c, d, e, f]
594
595	#bac, bce, ecf, dbe
596	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
597
598	boundary = { (0, 0): 'First',
599	(0, 2): 'Second',
600	(2, 0): 'Third',
601	(2, 1): 'Fourth',
602	#(3, 1): 'Fifth', #Skip this
603	(3, 2): 'Sixth'}
604
605
606	mesh = Mesh(points, vertices) #, boundary)
607	mesh.check_integrity()
608
609	from anuga.config import default_boundary_tag
610	assert mesh.boundary[ (0, 0) ] == default_boundary_tag
611	assert mesh.boundary[ (0, 2) ] == default_boundary_tag
612	assert mesh.boundary[ (2, 0) ] == default_boundary_tag
613	assert mesh.boundary[ (2, 1) ] == default_boundary_tag
614	assert mesh.boundary[ (3, 1) ] == default_boundary_tag
615	assert mesh.boundary[ (3, 2) ] == default_boundary_tag
616
617
618	#Check enumeration
619	#for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
620	# b = -k-1
621	# assert mesh.neighbours[vol_id, edge_id] == b
622
623
624
625
626
627
628	def test_inputs(self):
629	a = [0.0, 0.0]
630	b = [0.0, 2.0]
631	c = [2.0,0.0]
632	d = [0.0, 4.0]
633	e = [2.0, 2.0]
634	f = [4.0,0.0]
635
636	points = [a, b, c, d, e, f]
637
638	#bac, bce, ecf, dbe
639	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
640
641	#Too few points
642	try:
643	mesh = Mesh([points[0]], vertices)
644	except AssertionError:
645	pass
646	else:
647	raise 'Should have raised an exception'
648
649	#Too few points - 1 element
650	try:
651	mesh = Mesh([points[0]], [vertices[0]])
652	except AssertionError:
653	pass
654	else:
655	raise 'Should have raised an exception'
656
657	#Wrong dimension of vertices
658	try:
659	mesh = Mesh(points, vertices[0])
660	except AssertionError:
661	pass
662	else:
663	raise 'Should have raised an exception'
664
665	#Unsubscriptable coordinates object raises exception
666	try:
667	mesh = Mesh(points[0], [vertices[0]])
668	except AssertionError:
669	pass
670	else:
671	raise 'Should have raised an exception'
672
673	#FIXME: This has been commented out pending a decision
674	#whether to allow partial boundary tags or not
675	#
676	#Not specifying all boundary tags
677	#try:
678	# mesh = Mesh(points, vertices, {(3,0): 'x'})
679	#except AssertionError:
680	# pass
681	#else:
682	# raise 'Should have raised an exception'
683
684	#Specifying wrong non existing segment
685	try:
686	mesh = Mesh(points, vertices, {(5,0): 'x'})
687	except AssertionError:
688	pass
689	else:
690	raise 'Should have raised an exception'
691
692
693
694
695	def test_internal_boundaries(self):
696	"""
697	get values based on triangle lists.
698	"""
699	from mesh_factory import rectangular
700	from Numeric import zeros, Float
701
702	#Create basic mesh
703	points, vertices, boundary = rectangular(1, 3)
704
705	# Add an internal boundary
706	boundary[(2,0)] = 'internal'
707	boundary[(1,0)] = 'internal'
708
709	#Create shallow water domain
710	domain = Mesh(points, vertices, boundary)
711	domain.build_tagged_elements_dictionary({'bottom':[0,1],
712	'top':[4,5],
713	'all':[0,1,2,3,4,5]})
714
715
716	def test_boundary_polygon(self):
717	from mesh_factory import rectangular
718	#from mesh import Mesh
719	from Numeric import zeros, Float
720
721	#Create basic mesh
722	points, vertices, boundary = rectangular(2, 2)
723	mesh = Mesh(points, vertices, boundary)
724
725
726	P = mesh.get_boundary_polygon()
727
728	assert len(P) == 8
729	assert allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],
730	[1.0, 0.5], [1.0, 1.0], [0.5, 1.0],
731	[0.0, 1.0], [0.0, 0.5]])
732	for p in points:
733	#print p, P
734	assert is_inside_polygon(p, P)
735
736
737	def test_boundary_polygon_II(self):
738	from Numeric import zeros, Float
739
740
741	#Points
742	a = [0.0, 0.0] #0
743	b = [0.0, 0.5] #1
744	c = [0.0, 1.0] #2
745	d = [0.5, 0.0] #3
746	e = [0.5, 0.5] #4
747	f = [1.0, 0.0] #5
748	g = [1.0, 0.5] #6
749	h = [1.0, 1.0] #7
750	i = [1.5, 0.5] #8
751
752	points = [a, b, c, d, e, f, g, h, i]
753
754	#dea, bae, bec, fgd,
755	#edg, ghe, gfi, gih
756	vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
757	[4,3,6], [6,7,4], [6,5,8], [6,8,7]]
758
759	mesh = Mesh(points, vertices)
760
761	mesh.check_integrity()
762
763	P = mesh.get_boundary_polygon()
764
765	assert len(P) == 8
766	assert allclose(P, [a, d, f, i, h, e, c, b])
767
768	for p in points:
769	#print p, P
770	assert is_inside_polygon(p, P)
771
772
773	def test_boundary_polygon_III(self):
774	"""Same as II but vertices ordered differently
775	"""
776
777	from Numeric import zeros, Float
778
779
780	#Points
781	a = [0.0, 0.0] #0
782	b = [0.0, 0.5] #1
783	c = [0.0, 1.0] #2
784	d = [0.5, 0.0] #3
785	e = [0.5, 0.5] #4
786	f = [1.0, 0.0] #5
787	g = [1.0, 0.5] #6
788	h = [1.0, 1.0] #7
789	i = [1.5, 0.5] #8
790
791	points = [a, b, c, d, e, f, g, h, i]
792
793	#edg, ghe, gfi, gih
794	#dea, bae, bec, fgd,
795	vertices = [[4,3,6], [6,7,4], [6,5,8], [6,8,7],
796	[3,4,0], [1,0,4], [1,4,2], [5,6,3]]
797
798
799	mesh = Mesh(points, vertices)
800	mesh.check_integrity()
801
802
803	P = mesh.get_boundary_polygon()
804
805	assert len(P) == 8
806	assert allclose(P, [a, d, f, i, h, e, c, b])
807
808	for p in points:
809	assert is_inside_polygon(p, P)
810
811
812	def test_boundary_polygon_IIIa(self):
813	"""test_boundary_polygon_IIIa - Check pathological situation where
814	one triangle has no neighbours. This may be the case if a mesh
815	is partitioned using pymetis.
816	"""
817
818	from Numeric import zeros, Float
819
820
821	#Points
822	a = [0.0, 0.0] #0
823	b = [0.0, 0.5] #1
824	c = [0.0, 1.0] #2
825	d = [0.5, 0.0] #3
826	e = [0.5, 0.5] #4
827	f = [1.0, 0.0] #5
828	g = [1.0, 0.5] #6
829	h = [1.0, 1.0] #7
830
831	# Add pathological triangle with no neighbours to an otherwise
832	# trivial mesh
833
834	points = [a, b, c, d, e, f, g, h]
835
836	#cbe, aeb, dea, fed, ghe (pathological triangle)
837	vertices = [[2,1,4], [0,4,1], [3,4,0], [5,4,3],
838	[6,7,4]]
839
840	mesh = Mesh(points, vertices)
841	mesh.check_integrity()
842
843	P = mesh.get_boundary_polygon(verbose=False)
844
845
846	assert len(P) == 9
847
848	# Note that point e appears twice!
849	assert allclose(P, [a, d, f, e, g, h, e, c, b])
850
851	for p in points:
852	msg = 'Point %s is not inside polygon %s'\
853	%(p, P)
854	assert is_inside_polygon(p, P), msg
855
856
857
858
859
860
861	def test_boundary_polygon_IV(self):
862	"""Reproduce test test_spatio_temporal_file_function_time
863	from test_util.py that looked as if it produced the wrong boundary
864	"""
865
866	from Numeric import zeros, Float
867	from mesh_factory import rectangular
868
869	#Create a domain to hold test grid
870	#(0:15, -20:10)
871	points, vertices, boundary =\
872	rectangular(4, 4, 15, 30, origin = (0, -20))
873
874	#####
875	mesh = Mesh(points, vertices)
876	mesh.check_integrity()
877
878	P = mesh.get_boundary_polygon()
879
880	#print P
881	assert len(P) == 16
882	for p in points:
883	assert is_inside_polygon(p, P)
884
885
886
887	#####
888	mesh = Mesh(points, vertices, boundary)
889	mesh.check_integrity()
890
891	P = mesh.get_boundary_polygon()
892
893
894	#print P, len(P)
895	assert len(P) == 16
896
897	for p in points:
898	assert is_inside_polygon(p, P)
899
900	#print mesh.statistics()
901
902
903
904	def test_boundary_polygon_V(self):
905	"""Create a discontinuous mesh (duplicate vertices)
906	and check that boundary is as expected
907
908	"""
909	from Numeric import zeros, Float
910
911
912	#Points
913	a = [0.0, 0.0] #0
914	b = [0.0, 0.5] #1
915	c = [0.0, 1.0] #2
916	d = [0.5, 0.0] #3
917	e = [0.5, 0.5] #4
918	f = [1.0, 0.0] #5
919	g = [1.0, 0.5] #6
920	h = [1.0, 1.0] #7
921	i = [1.5, 0.5] #8
922
923	#Duplicate points for triangles edg [4,3,6] (central) and
924	#gid [6,8,7] (top right boundary) to them disconnected
925	#from the others
926
927	e0 = [0.5, 0.5] #9
928	d0 = [0.5, 0.0] #10
929	g0 = [1.0, 0.5] #11
930	i0 = [1.5, 0.5] #12
931
932
933	points = [a, b, c, d, e, f, g, h, i, e0, d0, g0, i0]
934
935
936
937	#dea, bae, bec, fgd,
938	#edg, ghe, gfi, gih
939	#vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
940	# [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
941
942
943	#dea, bae, bec, fgd,
944	#e0d0g0, ghe, gfi, g0i0h
945	vertices = [[3,4,0], [1,0,4], [1,4,2], [5,6,3],
946	[9,10,11], [6,7,4], [6,5,8], [11,12,7]]
947
948	mesh = Mesh(points, vertices)
949
950	mesh.check_integrity()
951
952	P = mesh.get_boundary_polygon()
953
954	#print P
955
956	assert len(P) == 8
957	assert allclose(P, [a, d, f, i, h, e, c, b])
958	assert allclose(P, [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.5, 0.5), (1.0, 1.0), (0.5, 0.5), (0.0, 1.0), (0.0, 0.5)])
959
960
961	for p in points:
962	#print p, P
963	assert is_inside_polygon(p, P)
964
965
966
967	def test_boundary_polygon_VI(self):
968	"""test_boundary_polygon_VI(self)
969
970	Create a discontinuous mesh (duplicate vertices) from a real situation that failed
971	and check that boundary is as expected
972	"""
973
974
975	from anuga.utilities.polygon import plot_polygons
976
977	# First do the continuous version of mesh
978
979	points = [[ 6626.85400391, 0. ],
980	[ 0. , 38246.4140625 ],
981	[ 9656.2734375 , 68351.265625 ],
982	[ 20827.25585938, 77818.203125 ],
983	[ 32755.59375 , 58126.9765625 ],
984	[ 35406.3359375 , 79332.9140625 ],
985	[ 31998.23828125, 88799.84375 ],
986	[ 23288.65820313, 104704.296875 ],
987	[ 32187.57617188, 109816.4375 ],
988	[ 50364.08984375, 110763.1328125 ],
989	[ 80468.9453125 , 96184.0546875 ],
990	[ 86149.1015625 , 129886.34375 ],
991	[ 118715.359375 , 129886.34375 ],
992	[ 117768.6640625 , 85770.4296875 ],
993	[ 101485.5390625 , 45251.9453125 ],
994	[ 49985.4140625 , 2272.06396484],
995	[ 51737.94140625, 90559.2109375 ],
996	[ 56659.0703125 , 65907.6796875 ],
997	[ 75735.4765625 , 23762.00585938],
998	[ 52341.70703125, 38563.39453125]]
999
1000	##points = ensure_numeric(points, Int)/1000 # Simplify for ease of interpretation
1001
1002	triangles = [[19, 0,15],
1003	[ 2, 4, 3],
1004	[ 4, 2, 1],
1005	[ 1,19, 4],
1006	[15,18,19],
1007	[18,14,17],
1008	[19, 1, 0],
1009	[ 6, 8, 7],
1010	[ 8, 6,16],
1011	[10, 9,16],
1012	[17, 5, 4],
1013	[16,17,10],
1014	[17,19,18],
1015	[ 5,17,16],
1016	[10,14,13],
1017	[10,17,14],
1018	[ 8,16, 9],
1019	[12,11,10],
1020	[10,13,12],
1021	[19,17, 4],
1022	[16, 6, 5]]
1023
1024	mesh = Mesh(points, triangles)
1025	mesh.check_integrity()
1026	Pref = mesh.get_boundary_polygon()
1027
1028	#plot_polygons([ensure_numeric(Pref)], 'goodP')
1029
1030	for p in points:
1031	assert is_inside_polygon(p, Pref)
1032
1033
1034	# Then do the discontinuous version
1035	import warnings
1036	warnings.filterwarnings('ignore')
1037
1038
1039	points = [[ 52341.70703125, 38563.39453125],
1040	[ 6626.85400391, 0. ],
1041	[ 49985.4140625 , 2272.06396484],
1042	[ 9656.2734375 , 68351.265625 ],
1043	[ 32755.59375 , 58126.9765625 ],
1044	[ 20827.25585938, 77818.203125 ],
1045	[ 32755.59375 , 58126.9765625 ],
1046	[ 9656.2734375 , 68351.265625 ],
1047	[ 0. , 38246.4140625 ],
1048	[ 0. , 38246.4140625 ],
1049	[ 52341.70703125, 38563.39453125],
1050	[ 32755.59375 , 58126.9765625 ],
1051	[ 49985.4140625 , 2272.06396484],
1052	[ 75735.4765625 , 23762.00585938],
1053	[ 52341.70703125, 38563.39453125],
1054	[ 75735.4765625 , 23762.00585938],
1055	[ 101485.5390625 , 45251.9453125 ],
1056	[ 56659.0703125 , 65907.6796875 ],
1057	[ 52341.70703125, 38563.39453125],
1058	[ 0. , 38246.4140625 ],
1059	[ 6626.85400391, 0. ],
1060	[ 31998.23828125, 88799.84375 ],
1061	[ 32187.57617188, 109816.4375 ],
1062	[ 23288.65820313, 104704.296875 ],
1063	[ 32187.57617188, 109816.4375 ],
1064	[ 31998.23828125, 88799.84375 ],
1065	[ 51737.94140625, 90559.2109375 ],
1066	[ 80468.9453125 , 96184.0546875 ],
1067	[ 50364.08984375, 110763.1328125 ],
1068	[ 51737.94140625, 90559.2109375 ],
1069	[ 56659.0703125 , 65907.6796875 ],
1070	[ 35406.3359375 , 79332.9140625 ],
1071	[ 32755.59375 , 58126.9765625 ],
1072	[ 51737.94140625, 90559.2109375 ],
1073	[ 56659.0703125 , 65907.6796875 ],
1074	[ 80468.9453125 , 96184.0546875 ],
1075	[ 56659.0703125 , 65907.6796875 ],
1076	[ 52341.70703125, 38563.39453125],
1077	[ 75735.4765625 , 23762.00585938],
1078	[ 35406.3359375 , 79332.9140625 ],
1079	[ 56659.0703125 , 65907.6796875 ],
1080	[ 51737.94140625, 90559.2109375 ],
1081	[ 80468.9453125 , 96184.0546875 ],
1082	[ 101485.5390625 , 45251.9453125 ],
1083	[ 117768.6640625 , 85770.4296875 ],
1084	[ 80468.9453125 , 96184.0546875 ],
1085	[ 56659.0703125 , 65907.6796875 ],
1086	[ 101485.5390625 , 45251.9453125 ],
1087	[ 32187.57617188, 109816.4375 ],
1088	[ 51737.94140625, 90559.2109375 ],
1089	[ 50364.08984375, 110763.1328125 ],
1090	[ 118715.359375 , 129886.34375 ],
1091	[ 86149.1015625 , 129886.34375 ],
1092	[ 80468.9453125 , 96184.0546875 ],
1093	[ 80468.9453125 , 96184.0546875 ],
1094	[ 117768.6640625 , 85770.4296875 ],
1095	[ 118715.359375 , 129886.34375 ],
1096	[ 52341.70703125, 38563.39453125],
1097	[ 56659.0703125 , 65907.6796875 ],
1098	[ 32755.59375 , 58126.9765625 ],
1099	[ 51737.94140625, 90559.2109375 ],
1100	[ 31998.23828125, 88799.84375 ],
1101	[ 35406.3359375 , 79332.9140625 ]]
1102
1103	scaled_points = ensure_numeric(points, Int)/1000 # Simplify for ease of interpretation
1104
1105	triangles = [[ 0, 1, 2],
1106	[ 3, 4, 5],
1107	[ 6, 7, 8],
1108	[ 9,10,11],
1109	[12,13,14],
1110	[15,16,17],
1111	[18,19,20],
1112	[21,22,23],
1113	[24,25,26],
1114	[27,28,29],
1115	[30,31,32],
1116	[33,34,35],
1117	[36,37,38],
1118	[39,40,41],
1119	[42,43,44],
1120	[45,46,47],
1121	[48,49,50],
1122	[51,52,53],
1123	[54,55,56],
1124	[57,58,59],
1125	[60,61,62]]
1126
1127
1128	# First use scaled points for ease of debugging
1129	mesh = Mesh(scaled_points, triangles)
1130	mesh.check_integrity()
1131	P = mesh.get_boundary_polygon()
1132
1133	for p in scaled_points:
1134	assert is_inside_polygon(p, P)
1135
1136	# Then use original points and test
1137	mesh = Mesh(points, triangles)
1138	mesh.check_integrity()
1139	P = mesh.get_boundary_polygon()
1140
1141	for p in points:
1142	assert is_inside_polygon(p, P)
1143
1144	assert allclose(P, Pref)
1145
1146	def test_lone_vertices(self):
1147	a = [2.0, 1.0]
1148	b = [6.0, 2.0]
1149	c = [1.0, 3.0]
1150	d = [2.0, 4.0]
1151
1152	points = [a, b, c, d]
1153	vertices = [[0,1,2]]
1154
1155	mesh = Mesh(points, vertices)
1156	mesh.check_integrity()
1157	loners = mesh.get_lone_vertices()
1158	self.failUnless(loners==[3],
1159	'FAILED!')
1160
1161
1162	a = [2.0, 1.0]
1163	b = [6.0, 2.0]
1164	c = [1.0, 3.0]
1165	d = [2.0, 4.0]
1166
1167	points = [d, a, b, c]
1168	vertices = [[3,1,2]]
1169
1170	mesh = Mesh(points, vertices)
1171	mesh.check_integrity()
1172	loners = mesh.get_lone_vertices()
1173	self.failUnless(loners==[0],
1174	'FAILED!')
1175
1176	def test_mesh_get_boundary_polygon_with_georeferencing(self):
1177
1178	# test
1179	a = [0.0, 0.0]
1180	b = [4.0, 0.0]
1181	c = [0.0, 4.0]
1182
1183	absolute_points = [a, b, c]
1184	vertices = [[0,1,2]]
1185
1186	geo = Geo_reference(56,67,-56)
1187
1188	relative_points = geo.change_points_geo_ref(absolute_points)
1189
1190	#print 'Relative', relative_points
1191	#print 'Absolute', absolute_points
1192
1193	mesh = Mesh(relative_points, vertices, geo_reference=geo)
1194	boundary_polygon = mesh.get_boundary_polygon()
1195
1196	assert allclose(absolute_points, boundary_polygon)
1197
1198	def test_get_triangle_containing_point(self):
1199
1200	a = [0.0, 0.0]
1201	b = [0.0, 2.0]
1202	c = [2.0, 0.0]
1203	d = [0.0, 4.0]
1204	e = [2.0, 2.0]
1205	f = [4.0, 0.0]
1206
1207	points = [a, b, c, d, e, f]
1208	#bac, bce, ecf, dbe
1209	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1210	mesh = Mesh(points, vertices)
1211
1212	mesh.check_integrity()
1213
1214
1215	try:
1216	id = mesh.get_triangle_containing_point([3.0, 5.0])
1217	except:
1218	pass
1219	else:
1220	msg = 'Should have caught point outside polygon (Non)'
1221	raise Exception, msg
1222
1223	id = mesh.get_triangle_containing_point([0.5, 1.0])
1224	assert id == 0
1225
1226	id = mesh.get_triangle_containing_point([1.0, 3.0])
1227	assert id == 3
1228
1229	for i, point in enumerate(mesh.get_centroid_coordinates()):
1230	id = mesh.get_triangle_containing_point(point)
1231	assert id == i
1232
1233	def test_get_triangle_neighbours(self):
1234	a = [0.0, 0.0]
1235	b = [0.0, 2.0]
1236	c = [2.0,0.0]
1237	e = [2.0, 2.0]
1238	points = [a, b, c, e]
1239	vertices = [ [1,0,2], [1,2,3] ] #bac, bce
1240	mesh = Mesh(points, vertices)
1241
1242	neighbours = mesh.get_triangle_neighbours(0)
1243	assert allclose(neighbours, [-1,1,-1])
1244	neighbours = mesh.get_triangle_neighbours(-10)
1245	assert neighbours == []
1246	neighbours = mesh.get_triangle_neighbours(2)
1247	assert neighbours == []
1248
1249
1250	def test_get_intersecting_segments1(self):
1251	"""test_get_intersecting_segments(self):
1252
1253	Very simple test (horizontal lines)
1254	"""
1255
1256	# Build test mesh
1257
1258	# Create basic mesh
1259	# 9 points at (0,0), (0, 1), ..., (2,2)
1260	# 8 triangles enumerated from left bottom to right top.
1261	points, vertices, boundary = rectangular(2, 2, 2, 2)
1262	mesh = Mesh(points, vertices, boundary)
1263
1264	# Very simple horizontal line intersecting
1265	#
1266
1267
1268	for y_line in [0.1, 0.2, 0.314159, 0.41, 0.6, 0.99, 1.01, 1.5, 1.77, 1.9]:
1269	if y_line < 1:
1270	ceiling = 1
1271	floor = 0
1272	intersected_triangles = [0,1,4,5]
1273	elif y_line > 1:
1274	ceiling = 2
1275	floor = 1
1276	intersected_triangles = [2,3,6,7]
1277	else:
1278	raise Exception, 'this test is not for parallel lines'
1279
1280
1281	line = [[-1,y_line], [3,y_line]]
1282
1283	L = mesh.get_intersecting_segments(line)
1284	assert len(L) == 4
1285
1286
1287
1288	# Check all normals point straight down etc
1289	total_length = 0
1290	for x in L:
1291	if x.triangle_id % 2 == 0:
1292	assert allclose(x.length, ceiling-y_line)
1293	else:
1294	assert allclose(x.length, y_line-floor)
1295
1296
1297	assert allclose(x.normal, [0,-1])
1298
1299	assert allclose(x.segment[1][0], x.segment[0][0] + x.length)
1300	assert allclose(x.segment[0][1], y_line)
1301	assert allclose(x.segment[1][1], y_line)
1302
1303	assert x.triangle_id in intersected_triangles
1304
1305	total_length += x.length
1306
1307	msg = 'Segments do not add up'
1308	assert allclose(total_length, 2), msg
1309
1310
1311	def test_get_intersecting_segments_coinciding(self):
1312	"""test_get_intersecting_segments_coinciding(self):
1313
1314	Test that lines coinciding with triangle edges work.
1315	"""
1316
1317	# Build test mesh
1318
1319	# Create basic mesh
1320	# 9 points at (0,0), (0, 1), ..., (2,2)
1321	# 8 triangles enumerated from left bottom to right top.
1322	points, vertices, boundary = rectangular(2, 2, 2, 2)
1323	mesh = Mesh(points, vertices, boundary)
1324	intersected_triangles = [1,5]
1325
1326	# Very simple horizontal line intersecting
1327	#
1328
1329	y_line = 1.0
1330
1331	line = [[-1,y_line], [3,y_line]]
1332
1333	L = mesh.get_intersecting_segments(line)
1334
1335
1336	msg = 'Only two triangles should be returned'
1337	assert len(L) == 2, msg
1338
1339
1340	# Check all
1341	total_length = 0
1342	for x in L:
1343	assert allclose(x.length, 1.0)
1344	assert allclose(x.normal, [0,-1])
1345
1346	assert allclose(x.segment[1][0], x.segment[0][0] + x.length)
1347	assert allclose(x.segment[0][1], y_line)
1348	assert allclose(x.segment[1][1], y_line)
1349
1350
1351
1352	assert x.triangle_id in intersected_triangles
1353
1354	total_length += x.length
1355
1356	msg = 'Segments do not add up'
1357	assert allclose(total_length, 2), msg
1358
1359
1360	def test_get_intersecting_segments_partially_coinciding(self):
1361	"""test_get_intersecting_segments_partially_coinciding(self):
1362
1363	Test that line coinciding with triangle edges work.
1364	But this ones only coincide with parts of the edge.
1365	"""
1366
1367	# Build test mesh
1368
1369	# Create basic mesh
1370	# 9 points at (0,0), (0, 1), ..., (2,2)
1371	# 8 triangles enumerated from left bottom to right top.
1372	points, vertices, boundary = rectangular(2, 2, 2, 2)
1373	mesh = Mesh(points, vertices, boundary)
1374	intersected_triangles = [1,5]
1375
1376	# Horizontal line intersecting along center but stopping
1377	# parway through second triangle's edge
1378	#
1379
1380	y_line = 1.0
1381
1382	#line = [[0, y_line], [2, y_line]]
1383	line = [[0, y_line], [1.5, y_line]]
1384
1385	L = mesh.get_intersecting_segments(line)
1386	#for x in L:
1387	# print x
1388
1389	msg = 'Two triangles should be returned'
1390	assert len(L) == 2, msg
1391
1392
1393	# Check all
1394	total_length = 0
1395	for x in L:
1396	if x.triangle_id == 1:
1397	assert allclose(x.length, 1)
1398	assert allclose(x.normal, [0, -1])
1399
1400	if x.triangle_id == 5:
1401	assert allclose(x.length, 0.5)
1402	assert allclose(x.normal, [0, -1])
1403
1404
1405	assert x.triangle_id in intersected_triangles
1406
1407	total_length += x.length
1408
1409	msg = 'Segments do not add up'
1410	assert allclose(total_length, 1.5), msg
1411
1412
1413
1414	def test_get_intersecting_segments2(self):
1415	"""test_get_intersecting_segments(self):
1416
1417	Lines with a slope
1418	"""
1419
1420	s2 = sqrt(2.0)/2
1421
1422
1423	# Build test mesh
1424
1425	# Create basic mesh
1426	# 9 points at (0,0), (0, 1), ..., (2,2)
1427	# 8 triangles enumerated from left bottom to right top.
1428	points, vertices, boundary = rectangular(2, 2, 2, 2)
1429	mesh = Mesh(points, vertices, boundary)
1430
1431
1432	# Diagonal cutting through a vertex and hypothenuses
1433	line = [[0, 2], [2, 0]]
1434	intersected_triangles = [3,2,5,4]
1435
1436	L = mesh.get_intersecting_segments(line)
1437	assert len(L) == 4
1438
1439	#print L
1440
1441	# Check all segments
1442	total_length = 0
1443	for i, x in enumerate(L):
1444	assert allclose(x.length, s2)
1445	assert allclose(x.normal, [-s2, -s2])
1446	assert allclose(sum(x.normal**2), 1)
1447
1448	assert x.triangle_id in intersected_triangles
1449
1450	total_length += x.length
1451
1452	msg = 'Segments do not add up'
1453	assert allclose(total_length, 4*s2), msg
1454
1455
1456	# Diagonal cutting through a vertex and hypothenuses (reversed)
1457	line = [[2, 0], [0, 2]]
1458	intersected_triangles = [3,2,5,4]
1459
1460	L = mesh.get_intersecting_segments(line)
1461	assert len(L) == 4
1462
1463	#print L
1464
1465	# Check all segments
1466	total_length = 0
1467	for i, x in enumerate(L):
1468	assert allclose(x.length, s2)
1469	assert allclose(x.normal, [s2, s2])
1470	assert allclose(sum(x.normal**2), 1)
1471
1472	assert x.triangle_id in intersected_triangles
1473
1474	total_length += x.length
1475
1476	msg = 'Segments do not add up'
1477	assert allclose(total_length, 4*s2), msg
1478
1479
1480
1481	# Diagonal coinciding with hypothenuses
1482	line = [[2, 2], [0, 0]]
1483	intersected_triangles = [6,0]
1484
1485	L = mesh.get_intersecting_segments(line)
1486	assert len(L) == 2
1487
1488	#print L
1489
1490	# Check all segments
1491	total_length = 0
1492	for i, x in enumerate(L):
1493	assert allclose(x.length, 2*s2)
1494	assert allclose(x.normal, [-s2, s2])
1495	assert allclose(sum(x.normal**2), 1)
1496
1497	assert x.triangle_id in intersected_triangles
1498
1499	total_length += x.length
1500
1501	msg = 'Segments do not add up'
1502	assert allclose(total_length, 4*s2), msg
1503
1504
1505	# Diagonal coinciding with hypothenuses (reversed)
1506	line = [[0, 0], [2, 2]]
1507	intersected_triangles = [6,0]
1508
1509	L = mesh.get_intersecting_segments(line)
1510	assert len(L) == 2
1511
1512	#print L
1513
1514	# Check all segments
1515	total_length = 0
1516	for i, x in enumerate(L):
1517	assert allclose(x.length, 2*s2)
1518	assert allclose(x.normal, [s2, -s2])
1519	assert allclose(sum(x.normal**2), 1)
1520
1521	assert x.triangle_id in intersected_triangles
1522
1523	total_length += x.length
1524
1525	msg = 'Segments do not add up'
1526	assert allclose(total_length, 4*s2), msg
1527
1528
1529
1530	# line with slope [1, -1] cutting through vertices of tri 7 and 6
1531	line = [[1, 2], [2, 1]]
1532	intersected_triangles = [7,6]
1533
1534	L = mesh.get_intersecting_segments(line)
1535	assert len(L) == 2
1536
1537	#print L
1538
1539	# Check all segments
1540	total_length = 0
1541	for i, x in enumerate(L):
1542	assert allclose(x.length, s2)
1543	assert allclose(x.normal, [-s2, -s2])
1544	assert allclose(sum(x.normal**2), 1)
1545
1546	assert x.triangle_id in intersected_triangles
1547
1548	total_length += x.length
1549
1550	msg = 'Segments do not add up'
1551	assert allclose(total_length, 2*s2), msg
1552
1553
1554	# Arbitrary line with slope [1, -1] cutting through tri 7 and 6
1555	line = [[1.1, 2], [2.1, 1]]
1556	intersected_triangles = [7,6]
1557
1558	L = mesh.get_intersecting_segments(line)
1559	assert len(L) == 2
1560
1561	# Check all segments
1562	total_length = 0
1563	for i, x in enumerate(L):
1564	assert allclose(x.normal, [-s2, -s2])
1565	assert allclose(sum(x.normal**2), 1)
1566
1567	msg = 'Triangle %d' %x.triangle_id + ' is not in %s' %(intersected_triangles)
1568	assert x.triangle_id in intersected_triangles, msg
1569
1570
1571
1572	def test_get_intersecting_segments3(self):
1573	"""test_get_intersecting_segments(self):
1574
1575	Check that line can stop inside a triangle
1576
1577	"""
1578
1579
1580
1581	s2 = sqrt(2.0)/2
1582
1583
1584	# Build test mesh
1585
1586	# Create basic mesh
1587	# 9 points at (0,0), (0, 1), ..., (2,2)
1588	# 8 triangles enumerated from left bottom to right top.
1589	points, vertices, boundary = rectangular(2, 2, 2, 2)
1590	mesh = Mesh(points, vertices, boundary)
1591
1592
1593	# Line cutting through one triangle and ending on its edge
1594	line = [[0.5, 3], [0.5, 1.5]]
1595	intersected_triangles = [3]
1596
1597	L = mesh.get_intersecting_segments(line)
1598	assert len(L) == 1
1599	assert L[0].triangle_id == 3
1600	assert allclose(L[0].length, 0.5)
1601	assert allclose(L[0].normal, [-1,0])
1602
1603
1604
1605	# Now try to shorten it so that its endpoint falls short of the far edge
1606	line = [[0.5, 3], [0.5, 1.6]]
1607	intersected_triangles = [3]
1608
1609	L = mesh.get_intersecting_segments(line)
1610	assert len(L) == 1
1611	assert L[0].triangle_id == 3
1612	assert allclose(L[0].length, 0.4)
1613	assert allclose(L[0].normal, [-1,0])
1614
1615	intersected_triangles = [3]
1616
1617	# Now the same, but with direction changed
1618	line = [[0.5, 3], [0.5, 1.6]]
1619	line = [[0.5, 1.6], [0.5, 3]]
1620	intersected_triangles = [3]
1621
1622	L = mesh.get_intersecting_segments(line)
1623	assert len(L) == 1
1624	assert L[0].triangle_id == 3
1625	assert allclose(L[0].length, 0.4)
1626	assert allclose(L[0].normal, [1,0])
1627
1628
1629
1630
1631	def test_get_intersecting_segments4(self):
1632	"""test_get_intersecting_segments(self):
1633
1634	Simple poly line
1635
1636	"""
1637
1638
1639
1640	s2 = sqrt(2.0)/2
1641
1642
1643	# Build test mesh
1644
1645	# Create basic mesh
1646	# 9 points at (0,0), (0, 1), ..., (2,2)
1647	# 8 triangles enumerated from left bottom to right top.
1648	points, vertices, boundary = rectangular(2, 2, 2, 2)
1649	mesh = Mesh(points, vertices, boundary)
1650
1651
1652	# Polyline with three segments cutting through mesh
1653	line = [[0.5, 3], [0.5, 1.5], [1,1]]
1654
1655	L = mesh.get_intersecting_segments(line)
1656	assert len(L) == 2
1657
1658	for x in L:
1659	if x.triangle_id == 3:
1660	assert allclose(x.length, 0.5)
1661	assert allclose(x.normal, [-1,0])
1662
1663	if x.triangle_id == 2:
1664	assert allclose(x.length, s2)
1665	assert allclose(x.normal, [-s2,-s2])
1666
1667
1668
1669	def test_get_intersecting_segments5(self):
1670	"""test_get_intersecting_segments(self):
1671
1672	More complex poly line
1673
1674	"""
1675
1676
1677
1678	s2 = sqrt(2.0)/2
1679
1680
1681	# Build test mesh
1682
1683	# Create basic mesh
1684	# 9 points at (0,0), (0, 1), ..., (2,2)
1685	# 8 triangles enumerated from left bottom to right top.
1686	points, vertices, boundary = rectangular(2, 2, 2, 2)
1687	mesh = Mesh(points, vertices, boundary)
1688
1689
1690	# Polyline with three segments cutting through mesh
1691	line = [[0.5, 3], [0.5, 1.5], [1.25, 0.75]]
1692
1693	L = mesh.get_intersecting_segments(line)
1694	assert len(L) == 3
1695
1696	for x in L:
1697	if x.triangle_id == 3:
1698	assert allclose(x.length, 0.5)
1699	assert allclose(x.normal, [-1,0])
1700
1701	if x.triangle_id == 2:
1702	msg = str(x.length)
1703	assert allclose(x.length, s2), msg
1704	assert allclose(x.normal, [-s2,-s2])
1705
1706	if x.triangle_id == 5:
1707	segvec = array([line[2][0]-1,
1708	line[2][1]-1])
1709	msg = str(x.length)
1710	assert allclose(x.length, sqrt(sum(segvec**2))), msg
1711	assert allclose(x.normal, [-s2,-s2])
1712
1713
1714	def test_get_intersecting_segments6(self):
1715	"""test_get_intersecting_segments(self):
1716
1717	Even more complex poly line, where line breaks within triangle 5
1718
1719	5 segments are returned even though only four triangles [3,2,5,6] are touched.
1720	Triangle 5 therefor has two segments in it.
1721
1722	"""
1723
1724
1725
1726	s2 = sqrt(2.0)/2
1727
1728
1729	# Build test mesh
1730
1731	# Create basic mesh
1732	# 9 points at (0,0), (0, 1), ..., (2,2)
1733	# 8 triangles enumerated from left bottom to right top.
1734	points, vertices, boundary = rectangular(2, 2, 2, 2)
1735	mesh = Mesh(points, vertices, boundary)
1736
1737
1738	# Polyline with three segments cutting through mesh
1739	line = [[0.5, 3], [0.5, 1.5], [1.25, 0.75], [2.25, 1.75]]
1740
1741	L = mesh.get_intersecting_segments(line)
1742	#for x in L:
1743	# print x
1744
1745	assert len(L) == 5
1746
1747	for x in L:
1748	if x.triangle_id == 3:
1749	assert allclose(x.length, 0.5)
1750	assert allclose(x.normal, [-1,0])
1751
1752	if x.triangle_id == 2:
1753	msg = str(x.length)
1754	assert allclose(x.length, s2), msg
1755	assert allclose(x.normal, [-s2,-s2])
1756
1757	if x.triangle_id == 5:
1758	if x.segment == ((1.0, 1.0), (1.25, 0.75)):
1759	segvec = array([line[2][0]-1,
1760	line[2][1]-1])
1761	msg = str(x.length)
1762	assert allclose(x.length, sqrt(sum(segvec**2))), msg
1763	assert allclose(x.normal, [-s2,-s2])
1764	elif x.segment == ((1.25, 0.75), (1.5, 1.0)):
1765	segvec = array([1.5-line[2][0],
1766	1.0-line[2][1]])
1767
1768	assert allclose(x.length, sqrt(sum(segvec**2))), msg
1769	assert allclose(x.normal, [s2,-s2])
1770	else:
1771	msg = 'Unknow segment: %s' %x.segment
1772	raise Exception, msg
1773
1774
1775
1776	if x.triangle_id == 6:
1777	assert allclose(x.normal, [s2,-s2])
1778	assert allclose(x.segment, ((1.5, 1.0), (2, 1.5)))
1779
1780
1781	# Internal test that sum of line segments add up
1782	# to length of input line
1783	#
1784	# Could be useful perhaps
1785	#
1786	#xi1 = line[1][0]
1787	#eta1 = line[1][1]
1788	#linevector = array([xi1-xi0, eta1-eta0])
1789	#linelength = sqrt(sum(linevector**2))
1790	#
1791	#segmentlength = 0
1792	#for segment in triangle_intersections:
1793	# vector = array([segment[1][0] - segment[0][0],
1794	# segment[1][1] - segment[0][1]])
1795	# length = sqrt(sum(vector**2))
1796	# segmentlength += length
1797	#
1798	#msg = 'Sum of intersecting segments do not add up'
1799	#assert allclose(segmentlength, linelength), msg
1800
1801
1802
1803
1804	def test_get_intersecting_segments7(self):
1805	"""test_get_intersecting_segments(self):
1806
1807	Check that line can stop inside a triangle - this is from
1808	flow throug a cross sections example in test_datamanager.
1809
1810	"""
1811
1812	# Build test mesh
1813	width = 5
1814	length = 100
1815	t_end = 1
1816	points, vertices, boundary = rectangular(length, width,
1817	length, width)
1818
1819	mesh = Mesh(points, vertices, boundary)
1820
1821
1822	# A range of partial lines
1823	x = length/2.
1824	for i in range(10):
1825	start_point = [length/2., i*width/10.]
1826	#print
1827	#print start_point
1828
1829	line = [start_point, [length/2., width]]
1830
1831	L = mesh.get_intersecting_segments(line)
1832
1833	if start_point[1] < 1:
1834	assert len(L) == 5
1835
1836
1837	total_length = 0
1838	for x in L:
1839	total_length += x.length
1840
1841
1842	ref_length = line[1][1] - line[0][1]
1843	#print ref_length, total_length
1844	assert allclose(total_length, ref_length)
1845
1846
1847
1848
1849	#-------------------------------------------------------------
1850	if __name__ == "__main__":
1851	#suite = unittest.makeSuite(Test_Mesh,'test_two_triangles')
1852	#suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_partially_coinciding')
1853	#suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments7')
1854	suite = unittest.makeSuite(Test_Mesh,'test')
1855	runner = unittest.TextTestRunner()
1856	runner.run(suite)
1857
1858
1859
1860

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

Download in other formats: