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

Last change on this file since 6145 was 6145, checked in by rwilson, 15 years ago

Change Numeric imports to general form - ready to change to NumPy?.

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