Context Navigation

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

Last change on this file since 3945 was 3945, checked in by ole, 17 years ago
One large step towards major cleanup. This has mainly to do with the way vertex coordinates are handled internally.
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: