Context Navigation

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

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

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

Download in other formats: