source: inundation/pyvolution/mesh_factory.py @ 1882

Last change on this file since 1882 was 1762, checked in by steve, 19 years ago
File size: 17.9 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 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 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 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 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 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 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 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 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
Note: See TracBrowser for help on using the repository browser.