source: anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/mesh_factory.py @ 5899

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

Initial NumPy? changes (again!).

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