source: anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py @ 3560

Last change on this file since 3560 was 3560, checked in by ole, 18 years ago

Renamed pyvolution to abstract_2d_finite_volumes. This is
one step towards pulling pyvolution apart. More to follow.
All unit tests pass and most examples fixed up.




File size: 20.3 KB
Line 
1"""Library of standard meshes and facilities for reading various
2mesh file formats
3"""
4
5
6def rectangular_old(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
7
8    """Setup a rectangular grid of triangles
9    with m+1 by n+1 grid points
10    and side lengths len1, len2. If side lengths are omitted
11    the mesh defaults to the unit square.
12
13    len1: x direction (left to right)
14    len2: y direction (bottom to top)
15
16    Return to lists: points and elements suitable for creating a Mesh or
17    FVMesh object, e.g. Mesh(points, elements)
18    """
19
20    from anuga.config import epsilon
21
22    deltax = float(len1)/m
23    deltay = float(len2)/n
24
25    #Dictionary of vertex objects
26    vertices = {}
27    points = []
28
29    for i in range(m+1):
30        for j in range(n+1):
31            vertices[i,j] = len(points)
32            points.append([i*delta1 + origin[0], j*delta2 + origin[1]])
33
34
35    #Construct 2 triangles per rectangular element and assign tags to boundary
36    elements = []
37    boundary = {}
38    for i in range(m):
39        for j in range(n):
40            v1 = vertices[i,j+1]
41            v2 = vertices[i,j]
42            v3 = vertices[i+1,j+1]
43            v4 = vertices[i+1,j]
44
45            #Update boundary dictionary and create elements
46            if i == m-1:
47                boundary[(len(elements), 2)] = 'right'
48            if j == 0:
49                boundary[(len(elements), 1)] = 'bottom'
50            elements.append([v4,v3,v2]) #Lower element
51
52            if i == 0:
53                boundary[(len(elements), 2)] = 'left'
54            if j == n-1:
55                boundary[(len(elements), 1)] = 'top'
56            elements.append([v1,v2,v3]) #Upper element
57
58    return points, elements, boundary
59
60def rectangular(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
61
62    """Setup a rectangular grid of triangles
63    with m+1 by n+1 grid points
64    and side lengths len1, len2. If side lengths are omitted
65    the mesh defaults to the unit square.
66
67    len1: x direction (left to right)
68    len2: y direction (bottom to top)
69
70    Return to lists: points and elements suitable for creating a Mesh or
71    FVMesh object, e.g. Mesh(points, elements)
72    """
73
74    from anuga.config import epsilon
75    from Numeric import zeros, Float, Int
76
77    delta1 = float(len1)/m
78    delta2 = float(len2)/n
79
80    #Calculate number of points
81    Np = (m+1)*(n+1)
82
83    class Index:
84
85        def __init__(self, n,m):
86            self.n = n
87            self.m = m
88
89        def __call__(self, i,j):
90            return j+i*(self.n+1)
91
92
93    index = Index(n,m)
94
95    points = zeros( (Np,2), Float)
96
97    for i in range(m+1):
98        for j in range(n+1):
99
100            points[index(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
101
102    #Construct 2 triangles per rectangular element and assign tags to boundary
103    #Calculate number of triangles
104    Nt = 2*m*n
105
106
107    elements = zeros( (Nt,3), Int)
108    boundary = {}
109    nt = -1
110    for i in range(m):
111        for j in range(n):
112            nt = nt + 1
113            i1 = index(i,j+1)
114            i2 = index(i,j)
115            i3 = index(i+1,j+1)
116            i4 = index(i+1,j)
117
118
119            #Update boundary dictionary and create elements
120            if i == m-1:
121                boundary[nt, 2] = 'right'
122            if j == 0:
123                boundary[nt, 1] = 'bottom'
124            elements[nt,:] = [i4,i3,i2] #Lower element
125            nt = nt + 1
126
127            if i == 0:
128                boundary[nt, 2] = 'left'
129            if j == n-1:
130                boundary[nt, 1] = 'top'
131            elements[nt,:] = [i1,i2,i3] #Upper element
132
133    return points, elements, boundary
134
135
136def rectangular_cross(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
137
138    """Setup a rectangular grid of triangles
139    with m+1 by n+1 grid points
140    and side lengths len1, len2. If side lengths are omitted
141    the mesh defaults to the unit square.
142
143    len1: x direction (left to right)
144    len2: y direction (bottom to top)
145
146    Return to lists: points and elements suitable for creating a Mesh or
147    FVMesh object, e.g. Mesh(points, elements)
148    """
149
150    from anuga.config import epsilon
151    from Numeric import zeros, Float, Int
152
153    delta1 = float(len1)/m
154    delta2 = float(len2)/n
155
156    #Dictionary of vertex objects
157    vertices = {}
158    points = []
159
160    for i in range(m+1):
161        for j in range(n+1):
162            vertices[i,j] = len(points)
163            points.append([delta1*i + origin[0], delta2*+ origin[1]])
164
165    # Construct 4 triangles per element
166    elements = []
167    boundary = {}
168    for i in range(m):
169        for j in range(n):
170            v1 = vertices[i,j+1]
171            v2 = vertices[i,j]
172            v3 = vertices[i+1,j+1]
173            v4 = vertices[i+1,j]
174            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
175            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
176
177            # Create centre point
178            v5 = len(points)
179            points.append([x, y])
180
181            #Create left triangle
182            if i == 0:
183                boundary[(len(elements), 1)] = 'left'
184            elements.append([v2,v5,v1])
185
186            #Create bottom triangle
187            if j == 0:
188                boundary[(len(elements), 1)] = 'bottom'
189            elements.append([v4,v5,v2])
190
191            #Create right triangle
192            if i == m-1:
193                boundary[(len(elements), 1)] = 'right'
194            elements.append([v3,v5,v4])
195
196            #Create top triangle
197            if j == n-1:
198                boundary[(len(elements), 1)] = 'top'
199            elements.append([v1,v5,v3])
200
201
202    return points, elements, boundary
203
204def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
205    """Setup a oblique grid of triangles
206    with m segments in the x-direction and n segments in the y-direction
207
208    """
209
210    from Numeric import array
211    import math
212
213    from anuga.config import epsilon
214
215
216    deltax = lenx/float(m)
217    deltay = leny/float(n)
218    a = 0.75*lenx*math.tan(theta/180.*math.pi)
219    x1 = lenx
220    y1 = 0
221    x2 = lenx
222    y2 = leny
223    x3 = 0.25*lenx
224    y3 = leny
225    x4 = x3
226    y4 = 0
227    a2 = a/(x1-x4)
228    a1 = -a2*x4
229    a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3)
230    a3 = 1. - (a1 + a2*x3)/y3 - a4*x3
231
232    # Dictionary of vertex objects
233    vertices = {}
234    points = []
235
236    for i in range(m+1):
237        x = deltax*i
238        for j in range(n+1):
239            y = deltay*j
240            if x > 0.25*lenx:
241                y = a1 + a2*x + a3*y + a4*x*y
242
243            vertices[i,j] = len(points)
244            points.append([x + origin[0], y + origin[1]])
245
246    # Construct 2 triangles per element
247    elements = []
248    boundary = {}
249    for i in range(m):
250        for j in range(n):
251            v1 = vertices[i,j+1]
252            v2 = vertices[i,j]
253            v3 = vertices[i+1,j+1]
254            v4 = vertices[i+1,j]
255
256            #Update boundary dictionary and create elements
257            if i == m-1:
258                boundary[(len(elements), 2)] = 'right'
259            if j == 0:
260                boundary[(len(elements), 1)] = 'bottom'
261            elements.append([v4,v3,v2]) #Lower
262
263            if i == 0:
264                boundary[(len(elements), 2)] = 'left'
265            if j == n-1:
266                boundary[(len(elements), 1)] = 'top'
267
268            elements.append([v1,v2,v3]) #Upper
269
270    return points, elements, boundary
271
272
273def circular(m, n, radius=1.0, center = (0.0, 0.0)):
274    """Setup a circular grid of triangles with m concentric circles and
275    with n radial segments. If radius is are omitted the mesh defaults to
276    the unit circle radius.
277
278    radius: radius of circle
279
280    #FIXME: The triangles become degenerate for large values of m or n.
281    """
282
283
284
285    from math import pi, cos, sin
286
287    radius = float(radius) #Ensure floating point format
288
289    #Dictionary of vertex objects and list of points
290    vertices = {}
291    points = [[0.0, 0.0]] #Center point
292    vertices[0, 0] = 0
293
294    for i in range(n):
295        theta = 2*i*pi/n
296        x = cos(theta)
297        y = sin(theta)
298        for j in range(1,m+1):
299            delta = j*radius/m
300            vertices[i,j] = len(points)
301            points.append([delta*x, delta*y])
302
303    #Construct 2 triangles per element
304    elements = []
305    for i in range(n):
306        for j in range(1,m):
307
308            i1 = (i + 1) % n  #Wrap around
309
310            v1 = vertices[i,j+1]
311            v2 = vertices[i,j]
312            v3 = vertices[i1,j+1]
313            v4 = vertices[i1,j]
314
315            elements.append([v4,v2,v3]) #Lower
316            elements.append([v1,v3,v2]) #Upper
317
318
319    #Do the center
320    v1 = vertices[0,0]
321    for i in range(n):
322        i1 = (i + 1) % n  #Wrap around
323        v2 = vertices[i,1]
324        v3 = vertices[i1,1]
325
326        elements.append([v1,v2,v3]) #center
327
328    return points, elements
329
330
331def from_polyfile(name):
332    """Read mesh from .poly file, an obj like file format
333    listing first vertex coordinates and then connectivity
334    """
335
336    from anuga.abstract_2d_finite_volumes.util import anglediff
337    from math import pi
338    import os.path
339    root, ext = os.path.splitext(name)
340
341    if ext == 'poly':
342        filename = name
343    else:
344        filename = name + '.poly'
345
346
347    fid = open(filename)
348
349    points = []    #x, y
350    values = []    #z
351    ##vertex_values = []    #Repeated z
352    triangles = [] #v0, v1, v2
353
354    lines = fid.readlines()
355
356    keyword = lines[0].strip()
357    msg = 'First line in .poly file must contain the keyword: POINTS'
358    assert keyword == 'POINTS', msg
359
360    offending = 0
361    i = 1
362    while keyword == 'POINTS':
363        line = lines[i].strip()
364        i += 1
365
366        if line == 'POLYS':
367            keyword = line
368            break
369
370        fields = line.split(':')
371        assert int(fields[0]) == i-1, 'Point indices not consecutive'
372
373        #Split the three floats
374        xyz = fields[1].split()
375
376        x = float(xyz[0])
377        y = float(xyz[1])
378        z = float(xyz[2])
379
380        points.append([x, y])
381        values.append(z)
382
383
384    k = i
385    while keyword == 'POLYS':
386        line = lines[i].strip()
387        i += 1
388
389        if line == 'END':
390            keyword = line
391            break
392
393
394        fields = line.split(':')
395        assert int(fields[0]) == i-k, 'Poly indices not consecutive'
396
397        #Split the three indices
398        vvv = fields[1].split()
399
400        i0 = int(vvv[0])-1
401        i1 = int(vvv[1])-1
402        i2 = int(vvv[2])-1
403
404        #Check for and exclude degenerate areas
405        x0 = points[i0][0]
406        y0 = points[i0][1]
407        x1 = points[i1][0]
408        y1 = points[i1][1]
409        x2 = points[i2][0]
410        y2 = points[i2][1]
411
412        area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
413        if area > 0:
414
415            #Ensure that points are arranged in counter clock-wise order
416            v0 = [x1-x0, y1-y0]
417            v1 = [x2-x1, y2-y1]
418            v2 = [x0-x2, y0-y2]
419
420            a0 = anglediff(v1, v0)
421            a1 = anglediff(v2, v1)
422            a2 = anglediff(v0, v2)
423
424
425            if a0 < pi and a1 < pi and a2 < pi:
426                #all is well
427                j0 = i0
428                j1 = i1
429                j2 = i2
430            else:
431                #Swap two vertices
432                j0 = i1
433                j1 = i0
434                j2 = i2
435
436            triangles.append([j0, j1, j2])
437            ##vertex_values.append([values[j0], values[j1], values[j2]])
438        else:
439            offending +=1
440
441    print 'Removed %d offending triangles out of %d' %(offending, len(lines))
442    return points, triangles, values
443
444
445
446def strang_mesh(filename):
447    """Read Strang generated mesh.
448    """
449
450    from math import pi
451    from anuga.abstract_2d_finite_volumes.util import anglediff
452
453
454    fid = open(filename)
455    points = []    # List of x, y coordinates
456    triangles = [] # List of vertex ids as listed in the file
457
458    for line in fid.readlines():
459        fields = line.split()
460        if len(fields) == 2:
461            # we are reading vertex coordinates
462            points.append([float(fields[0]), float(fields[1])])
463        elif len(fields) == 3:
464            # we are reading triangle point id's (format ae+b)
465            triangles.append([int(float(fields[0]))-1,
466                              int(float(fields[1]))-1,
467                              int(float(fields[2]))-1])
468        else:
469            raise 'wrong format in ' + filename
470
471    elements = [] #Final list of elements
472
473    for t in triangles:
474        #Get vertex coordinates
475        v0 = t[0]
476        v1 = t[1]
477        v2 = t[2]
478
479        x0 = points[v0][0]
480        y0 = points[v0][1]
481        x1 = points[v1][0]
482        y1 = points[v1][1]
483        x2 = points[v2][0]
484        y2 = points[v2][1]
485
486        #Check that points are arranged in counter clock-wise order
487        vec0 = [x1-x0, y1-y0]
488        vec1 = [x2-x1, y2-y1]
489        vec2 = [x0-x2, y0-y2]
490
491        a0 = anglediff(vec1, vec0)
492        a1 = anglediff(vec2, vec1)
493        a2 = anglediff(vec0, vec2)
494
495        if a0 < pi and a1 < pi and a2 < pi:
496            elements.append([v0, v1, v2])
497        else:
498            elements.append([v0, v2, v1])
499
500    return points, elements
501
502# #Map from edge number to indices of associated vertices
503# edge_map = ((1,2), (0,2), (0,1))
504
505def contracting_channel(m, n, W_upstream = 1., W_downstream = 0.75,
506                            L_1 = 5.0, L_2 = 2.0, L_3 = 10, origin = (0.0, 0.0)):
507    """Setup a contracting channel grid of triangles
508    with m segments in the x-direction and n segments in the y-direction
509
510    """
511
512    from Numeric import array
513    import math
514
515    from anuga.config import epsilon
516
517
518    lenx = L_1 + L_2 + L_3
519    leny = W_upstream
520    deltax = lenx/float(m)
521    deltay = leny/float(n)
522
523    x1 = 0
524    y1 = 0
525    x2 = L_1
526    y2 = 0
527    x3 = L_1 + L_2
528    y3 = (W_upstream - W_downstream)/2
529    x4 = L_1 + L_2 + L_3
530    y4 = y3
531    x5 = x4
532    y5 = y4 + W_downstream
533    x6 = L_1 + L_2
534    y6 = y5
535    x7 = L_1
536    y7 = W_upstream
537    x8 = 0
538    y8 = W_upstream
539    a1 = 0
540    a2 = (W_upstream - W_downstream)/(2*L_2)
541    a3 = 1
542    a4 = (W_downstream - W_upstream)/(L_2*W_upstream)
543
544    # Dictionary of vertex objects
545    vertices = {}
546    points = []
547
548    for i in range(m+1):
549        x = deltax*i
550        for j in range(n+1):
551            y = deltay*j
552            if x > L_1 and x <= (L_1 + L_2):
553                y = a1 + a2*(x - L_1) + a3*y + a4*(x - L_1)*y
554            elif x > L_1 + L_2:
555                y = (W_upstream - W_downstream)/2 + deltay*j*W_downstream/W_upstream
556
557            vertices[i,j] = len(points)
558            points.append([x + origin[0], y + origin[1]])
559
560    # Construct 2 triangles per element
561    elements = []
562    boundary = {}
563    for i in range(m):
564        for j in range(n):
565            v1 = vertices[i,j+1]
566            v2 = vertices[i,j]
567            v3 = vertices[i+1,j+1]
568            v4 = vertices[i+1,j]
569
570            #Update boundary dictionary and create elements
571            if i == m-1:
572                boundary[(len(elements), 2)] = 'right'
573            if j == 0:
574                boundary[(len(elements), 1)] = 'bottom'
575            elements.append([v4,v3,v2]) #Lower
576
577            if i == 0:
578                boundary[(len(elements), 2)] = 'left'
579            if j == n-1:
580                boundary[(len(elements), 1)] = 'top'
581
582            elements.append([v1,v2,v3]) #Upper
583
584    return points, elements, boundary
585
586
587def contracting_channel_cross(m, n, W_upstream = 1., W_downstream = 0.75,
588                              L_1 = 5.0, L_2 = 2.0, L_3 = 10, origin = (0.0, 0.0)):
589    """Setup a contracting channel grid of triangles
590    with m segments in the x-direction and n segments in the y-direction
591
592    """
593
594    from Numeric import array
595    import math
596
597    from anuga.config import epsilon
598
599
600    lenx = L_1 + L_2 + L_3
601    leny = W_upstream
602    deltax = lenx/float(m)
603    deltay = leny/float(n)
604
605    x1 = 0
606    y1 = 0
607    x2 = L_1
608    y2 = 0
609    x3 = L_1 + L_2
610    y3 = (W_upstream - W_downstream)/2
611    x4 = L_1 + L_2 + L_3
612    y4 = y3
613    x5 = x4
614    y5 = y4 + W_downstream
615    x6 = L_1 + L_2
616    y6 = y5
617    x7 = L_1
618    y7 = W_upstream
619    x8 = 0
620    y8 = W_upstream
621    a1 = 0
622    a2 = (W_upstream - W_downstream)/(2*L_2)
623    a3 = 1
624    a4 = (W_downstream - W_upstream)/(L_2*W_upstream)
625
626    # Dictionary of vertex objects
627    vertices = {}
628    points = []
629
630    for i in range(m+1):
631        x = deltax*i
632        for j in range(n+1):
633            y = deltay*j
634            if x > L_1 and x <= (L_1 + L_2):
635                y = a1 + a2*(x - L_1) + a3*y + a4*(x - L_1)*y
636            elif x > L_1 + L_2:
637                y = (W_upstream - W_downstream)/2 + deltay*j*W_downstream/W_upstream
638
639            vertices[i,j] = len(points)
640            points.append([x + origin[0], y + origin[1]])
641
642    # Construct 4 triangles per element
643    elements = []
644    boundary = {}
645    for i in range(m):
646        for j in range(n):
647            v1 = vertices[i,j+1]
648            v2 = vertices[i,j]
649            v3 = vertices[i+1,j+1]
650            v4 = vertices[i+1,j]
651            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
652            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
653            v5 = len(points)
654            points.append([x, y])
655
656            #Create left triangle
657            if i == 0:
658                boundary[(len(elements), 1)] = 'left'
659            elements.append([v2,v5,v1])
660
661            #Create bottom triangle
662            if j == 0:
663                boundary[(len(elements), 1)] = 'bottom'
664            elements.append([v4,v5,v2])
665
666            #Create right triangle
667            if i == m-1:
668                boundary[(len(elements), 1)] = 'right'
669            elements.append([v3,v5,v4])
670
671            #Create top triangle
672            if j == n-1:
673                boundary[(len(elements), 1)] = 'top'
674            elements.append([v1,v5,v3])
675
676
677    return points, elements, boundary
678
679
680
681
682def oblique_cross(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
683    """Setup a oblique grid of triangles
684    with m segments in the x-direction and n segments in the y-direction
685
686    """
687
688    from Numeric import array
689    import math
690
691    from anuga.config import epsilon
692
693
694    deltax = lenx/float(m)
695    deltay = leny/float(n)
696    a = 0.75*lenx*math.tan(theta/180.*math.pi)
697    x1 = lenx
698    y1 = 0
699    x2 = lenx
700    y2 = leny
701    x3 = 0.25*lenx
702    y3 = leny
703    x4 = x3
704    y4 = 0
705    a2 = a/(x1-x4)
706    a1 = -a2*x4
707    a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3)
708    a3 = 1. - (a1 + a2*x3)/y3 - a4*x3
709
710    # Dictionary of vertex objects
711    vertices = {}
712    points = []
713
714    for i in range(m+1):
715        x = deltax*i
716        for j in range(n+1):
717            y = deltay*j
718            if x > 0.25*lenx:
719                y = a1 + a2*x + a3*y + a4*x*y
720
721            vertices[i,j] = len(points)
722            points.append([x + origin[0], y + origin[1]])
723
724    # Construct 4 triangles per element
725    elements = []
726    boundary = {}
727    for i in range(m):
728        for j in range(n):
729            v1 = vertices[i,j+1]
730            v2 = vertices[i,j]
731            v3 = vertices[i+1,j+1]
732            v4 = vertices[i+1,j]
733            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
734            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
735            v5 = len(points)
736            points.append([x, y])
737
738            #Update boundary dictionary and create elements
739                        #Create left triangle
740            if i == 0:
741                boundary[(len(elements), 1)] = 'left'
742            elements.append([v2,v5,v1])
743
744            #Create bottom triangle
745            if j == 0:
746                boundary[(len(elements), 1)] = 'bottom'
747            elements.append([v4,v5,v2])
748
749            #Create right triangle
750            if i == m-1:
751                boundary[(len(elements), 1)] = 'right'
752            elements.append([v3,v5,v4])
753
754            #Create top triangle
755            if j == n-1:
756                boundary[(len(elements), 1)] = 'top'
757            elements.append([v1,v5,v3])
758
759
760    return points, elements, boundary
Note: See TracBrowser for help on using the repository browser.