Context Navigation

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

Last change on this file since 4851 was 4851, checked in by sexton, 17 years ago
minor update
File size: 36.2 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_points
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	#-------------------------------------------------------------
1234	if __name__ == "__main__":
1235	#suite = unittest.makeSuite(Test_Mesh,'test_get_triangle_containing_point')
1236	suite = unittest.makeSuite(Test_Mesh,'test')
1237	runner = unittest.TextTestRunner()
1238	runner.run(suite)
1239
1240
1241
1242

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

Download in other formats: