Context Navigation

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

Last change on this file since 4171 was 4171, checked in by ole, 17 years ago
Untabified a few culprits
File size: 34.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_IV(self):
813	"""Reproduce test test_spatio_temporal_file_function_time
814	from test_util.py that looked as if it produced the wrong boundary
815	"""
816
817	from Numeric import zeros, Float
818	from mesh_factory import rectangular
819
820	#Create a domain to hold test grid
821	#(0:15, -20:10)
822	points, vertices, boundary =\
823	rectangular(4, 4, 15, 30, origin = (0, -20))
824
825	#####
826	mesh = Mesh(points, vertices)
827	mesh.check_integrity()
828
829	P = mesh.get_boundary_polygon()
830
831	#print P
832	assert len(P) == 16
833	for p in points:
834	assert is_inside_polygon(p, P)
835
836
837
838	#####
839	mesh = Mesh(points, vertices, boundary)
840	mesh.check_integrity()
841
842	P = mesh.get_boundary_polygon()
843
844
845	#print P, len(P)
846	assert len(P) == 16
847
848	for p in points:
849	assert is_inside_polygon(p, P)
850
851	#print mesh.statistics()
852
853
854
855	def test_boundary_polygon_V(self):
856	"""Create a discontinuous mesh (duplicate vertices)
857	and check that boundary is as expected
858
859	"""
860	from Numeric import zeros, Float
861
862
863	#Points
864	a = [0.0, 0.0] #0
865	b = [0.0, 0.5] #1
866	c = [0.0, 1.0] #2
867	d = [0.5, 0.0] #3
868	e = [0.5, 0.5] #4
869	f = [1.0, 0.0] #5
870	g = [1.0, 0.5] #6
871	h = [1.0, 1.0] #7
872	i = [1.5, 0.5] #8
873
874	#Duplicate points for triangles edg [4,3,6] (central) and
875	#gid [6,8,7] (top right boundary) to them disconnected
876	#from the others
877
878	e0 = [0.5, 0.5] #9
879	d0 = [0.5, 0.0] #10
880	g0 = [1.0, 0.5] #11
881	i0 = [1.5, 0.5] #12
882
883
884	points = [a, b, c, d, e, f, g, h, i, e0, d0, g0, i0]
885
886
887
888	#dea, bae, bec, fgd,
889	#edg, ghe, gfi, gih
890	#vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
891	# [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
892
893
894	#dea, bae, bec, fgd,
895	#e0d0g0, ghe, gfi, g0i0h
896	vertices = [[3,4,0], [1,0,4], [1,4,2], [5,6,3],
897	[9,10,11], [6,7,4], [6,5,8], [11,12,7]]
898
899	mesh = Mesh(points, vertices)
900
901	mesh.check_integrity()
902
903	P = mesh.get_boundary_polygon()
904
905	#print P
906
907	assert len(P) == 8
908	assert allclose(P, [a, d, f, i, h, e, c, b])
909	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)])
910
911
912	for p in points:
913	#print p, P
914	assert is_inside_polygon(p, P)
915
916
917
918	def test_boundary_polygon_VI(self):
919	"""test_boundary_polygon_VI(self)
920
921	Create a discontinuous mesh (duplicate vertices) from a real situation that failed
922	and check that boundary is as expected
923	"""
924
925
926	from anuga.utilities.polygon import plot_polygons
927
928	# First do the continuous version of mesh
929
930	points = [[ 6626.85400391, 0. ],
931	[ 0. , 38246.4140625 ],
932	[ 9656.2734375 , 68351.265625 ],
933	[ 20827.25585938, 77818.203125 ],
934	[ 32755.59375 , 58126.9765625 ],
935	[ 35406.3359375 , 79332.9140625 ],
936	[ 31998.23828125, 88799.84375 ],
937	[ 23288.65820313, 104704.296875 ],
938	[ 32187.57617188, 109816.4375 ],
939	[ 50364.08984375, 110763.1328125 ],
940	[ 80468.9453125 , 96184.0546875 ],
941	[ 86149.1015625 , 129886.34375 ],
942	[ 118715.359375 , 129886.34375 ],
943	[ 117768.6640625 , 85770.4296875 ],
944	[ 101485.5390625 , 45251.9453125 ],
945	[ 49985.4140625 , 2272.06396484],
946	[ 51737.94140625, 90559.2109375 ],
947	[ 56659.0703125 , 65907.6796875 ],
948	[ 75735.4765625 , 23762.00585938],
949	[ 52341.70703125, 38563.39453125]]
950
951	##points = ensure_numeric(points, Int)/1000 # Simplify for ease of interpretation
952
953	triangles = [[19, 0,15],
954	[ 2, 4, 3],
955	[ 4, 2, 1],
956	[ 1,19, 4],
957	[15,18,19],
958	[18,14,17],
959	[19, 1, 0],
960	[ 6, 8, 7],
961	[ 8, 6,16],
962	[10, 9,16],
963	[17, 5, 4],
964	[16,17,10],
965	[17,19,18],
966	[ 5,17,16],
967	[10,14,13],
968	[10,17,14],
969	[ 8,16, 9],
970	[12,11,10],
971	[10,13,12],
972	[19,17, 4],
973	[16, 6, 5]]
974
975	mesh = Mesh(points, triangles)
976	mesh.check_integrity()
977	Pref = mesh.get_boundary_polygon()
978
979	#plot_polygons([ensure_numeric(Pref)], 'goodP')
980
981	for p in points:
982	assert is_inside_polygon(p, Pref)
983
984
985	# Then do the discontinuous version
986	import warnings
987	warnings.filterwarnings('ignore')
988
989
990	points = [[ 52341.70703125, 38563.39453125],
991	[ 6626.85400391, 0. ],
992	[ 49985.4140625 , 2272.06396484],
993	[ 9656.2734375 , 68351.265625 ],
994	[ 32755.59375 , 58126.9765625 ],
995	[ 20827.25585938, 77818.203125 ],
996	[ 32755.59375 , 58126.9765625 ],
997	[ 9656.2734375 , 68351.265625 ],
998	[ 0. , 38246.4140625 ],
999	[ 0. , 38246.4140625 ],
1000	[ 52341.70703125, 38563.39453125],
1001	[ 32755.59375 , 58126.9765625 ],
1002	[ 49985.4140625 , 2272.06396484],
1003	[ 75735.4765625 , 23762.00585938],
1004	[ 52341.70703125, 38563.39453125],
1005	[ 75735.4765625 , 23762.00585938],
1006	[ 101485.5390625 , 45251.9453125 ],
1007	[ 56659.0703125 , 65907.6796875 ],
1008	[ 52341.70703125, 38563.39453125],
1009	[ 0. , 38246.4140625 ],
1010	[ 6626.85400391, 0. ],
1011	[ 31998.23828125, 88799.84375 ],
1012	[ 32187.57617188, 109816.4375 ],
1013	[ 23288.65820313, 104704.296875 ],
1014	[ 32187.57617188, 109816.4375 ],
1015	[ 31998.23828125, 88799.84375 ],
1016	[ 51737.94140625, 90559.2109375 ],
1017	[ 80468.9453125 , 96184.0546875 ],
1018	[ 50364.08984375, 110763.1328125 ],
1019	[ 51737.94140625, 90559.2109375 ],
1020	[ 56659.0703125 , 65907.6796875 ],
1021	[ 35406.3359375 , 79332.9140625 ],
1022	[ 32755.59375 , 58126.9765625 ],
1023	[ 51737.94140625, 90559.2109375 ],
1024	[ 56659.0703125 , 65907.6796875 ],
1025	[ 80468.9453125 , 96184.0546875 ],
1026	[ 56659.0703125 , 65907.6796875 ],
1027	[ 52341.70703125, 38563.39453125],
1028	[ 75735.4765625 , 23762.00585938],
1029	[ 35406.3359375 , 79332.9140625 ],
1030	[ 56659.0703125 , 65907.6796875 ],
1031	[ 51737.94140625, 90559.2109375 ],
1032	[ 80468.9453125 , 96184.0546875 ],
1033	[ 101485.5390625 , 45251.9453125 ],
1034	[ 117768.6640625 , 85770.4296875 ],
1035	[ 80468.9453125 , 96184.0546875 ],
1036	[ 56659.0703125 , 65907.6796875 ],
1037	[ 101485.5390625 , 45251.9453125 ],
1038	[ 32187.57617188, 109816.4375 ],
1039	[ 51737.94140625, 90559.2109375 ],
1040	[ 50364.08984375, 110763.1328125 ],
1041	[ 118715.359375 , 129886.34375 ],
1042	[ 86149.1015625 , 129886.34375 ],
1043	[ 80468.9453125 , 96184.0546875 ],
1044	[ 80468.9453125 , 96184.0546875 ],
1045	[ 117768.6640625 , 85770.4296875 ],
1046	[ 118715.359375 , 129886.34375 ],
1047	[ 52341.70703125, 38563.39453125],
1048	[ 56659.0703125 , 65907.6796875 ],
1049	[ 32755.59375 , 58126.9765625 ],
1050	[ 51737.94140625, 90559.2109375 ],
1051	[ 31998.23828125, 88799.84375 ],
1052	[ 35406.3359375 , 79332.9140625 ]]
1053
1054	scaled_points = ensure_numeric(points, Int)/1000 # Simplify for ease of interpretation
1055
1056	triangles = [[ 0, 1, 2],
1057	[ 3, 4, 5],
1058	[ 6, 7, 8],
1059	[ 9,10,11],
1060	[12,13,14],
1061	[15,16,17],
1062	[18,19,20],
1063	[21,22,23],
1064	[24,25,26],
1065	[27,28,29],
1066	[30,31,32],
1067	[33,34,35],
1068	[36,37,38],
1069	[39,40,41],
1070	[42,43,44],
1071	[45,46,47],
1072	[48,49,50],
1073	[51,52,53],
1074	[54,55,56],
1075	[57,58,59],
1076	[60,61,62]]
1077
1078
1079	# First use scaled points for ease of debugging
1080	mesh = Mesh(scaled_points, triangles)
1081	mesh.check_integrity()
1082	P = mesh.get_boundary_polygon()
1083
1084	for p in scaled_points:
1085	assert is_inside_polygon(p, P)
1086
1087	# Then use original points and test
1088	mesh = Mesh(points, triangles)
1089	mesh.check_integrity()
1090	P = mesh.get_boundary_polygon()
1091
1092	for p in points:
1093	assert is_inside_polygon(p, P)
1094
1095	assert allclose(P, Pref)
1096
1097	def test_lone_vertices(self):
1098	a = [2.0, 1.0]
1099	b = [6.0, 2.0]
1100	c = [1.0, 3.0]
1101	d = [2.0, 4.0]
1102
1103	points = [a, b, c, d]
1104	vertices = [[0,1,2]]
1105
1106	mesh = Mesh(points, vertices)
1107	mesh.check_integrity()
1108	loners = mesh.get_lone_vertices()
1109	self.failUnless(loners==[3],
1110	'FAILED!')
1111
1112
1113	a = [2.0, 1.0]
1114	b = [6.0, 2.0]
1115	c = [1.0, 3.0]
1116	d = [2.0, 4.0]
1117
1118	points = [d, a, b, c]
1119	vertices = [[3,1,2]]
1120
1121	mesh = Mesh(points, vertices)
1122	mesh.check_integrity()
1123	loners = mesh.get_lone_vertices()
1124	self.failUnless(loners==[0],
1125	'FAILED!')
1126
1127	def test_mesh_get_boundary_polygon_with_georeferencing(self):
1128
1129	# test
1130	a = [0.0, 0.0]
1131	b = [4.0, 0.0]
1132	c = [0.0, 4.0]
1133
1134	absolute_points = [a, b, c]
1135	vertices = [[0,1,2]]
1136
1137	geo = Geo_reference(56,67,-56)
1138
1139	relative_points = geo.change_points_geo_ref(absolute_points)
1140
1141	#print 'Relative', relative_points
1142	#print 'Absolute', absolute_points
1143
1144	mesh = Mesh(relative_points, vertices, geo_reference=geo)
1145	boundary_polygon = mesh.get_boundary_polygon()
1146
1147	assert allclose(absolute_points, boundary_polygon)
1148
1149	#-------------------------------------------------------------
1150	if __name__ == "__main__":
1151	#suite = unittest.makeSuite(Test_Mesh,'test_mesh_get_boundary_polygon_with_georeferencing')
1152	suite = unittest.makeSuite(Test_Mesh,'test')
1153	runner = unittest.TextTestRunner()
1154	runner.run(suite)
1155
1156
1157
1158

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

Download in other formats: