source: inundation/pyvolution/mesh_factory_all.py @ 1869

Last change on this file since 1869 was 1531, checked in by chris, 19 years ago
File size: 13.5 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    #E = m*n*2        #Number of triangular elements
23    #P = (m+1)*(n+1)  #Number of initial vertices
24
25    delta1 = float(len1)/m
26    delta2 = float(len2)/n
27
28    #Dictionary of vertex objects
29    vertices = {}
30    points = []
31
32    for i in range(m+1):
33        for j in range(n+1):
34            vertices[i,j] = len(points)
35            points.append([i*delta1 + origin[0], j*delta2 + origin[1]])
36
37
38
39    #Construct 2 triangles per rectangular element and assign tags to boundary
40    elements = []
41    boundary = {}
42    for i in range(m):
43        for j in range(n):
44            v1 = vertices[i,j+1]
45            v2 = vertices[i,j]
46            v3 = vertices[i+1,j+1]
47            v4 = vertices[i+1,j]
48
49            #Update boundary dictionary and create elements
50            if i == m-1:
51                boundary[(len(elements), 2)] = 'right'
52            if j == 0:
53                boundary[(len(elements), 1)] = 'bottom'
54            elements.append([v4,v3,v2]) #Lower element
55
56            if i == 0:
57                boundary[(len(elements), 2)] = 'left'
58            if j == n-1:
59                boundary[(len(elements), 1)] = 'top'
60            elements.append([v1,v2,v3]) #Upper element
61
62    return points, elements, boundary
63
64def rectangular(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
65
66    """Setup a rectangular grid of triangles
67    with m+1 by n+1 grid points
68    and side lengths len1, len2. If side lengths are omitted
69    the mesh defaults to the unit square.
70
71    len1: x direction (left to right)
72    len2: y direction (bottom to top)
73
74    Return to lists: points and elements suitable for creating a Mesh or
75    FVMesh object, e.g. Mesh(points, elements)
76    """
77
78    from config import epsilon
79    from Numeric import zeros, Float, Int
80
81    delta1 = float(len1)/m
82    delta2 = float(len2)/n
83
84    #Calculate number of points
85    Np = (m+1)*(n+1)
86
87    class Index:
88
89        def __init__(self, n,m):
90            self.n = n
91            self.m = m
92
93        def __call__(self, i,j):
94            return j+i*(self.n+1)
95
96
97    index = Index(n,m)
98
99    points = zeros( (Np,2), Float)
100
101    for i in range(m+1):
102        for j in range(n+1):
103
104            points[index(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
105
106    #Construct 2 triangles per rectangular element and assign tags to boundary
107    #Calculate number of triangles
108    Nt = 2*m*n
109
110
111    elements = zeros( (Nt,3), Int)
112    boundary = {}
113    nt = -1
114    for i in range(m):
115        for j in range(n):
116            nt = nt + 1
117            i1 = index(i,j+1)
118            i2 = index(i,j)
119            i3 = index(i+1,j+1)
120            i4 = index(i+1,j)
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
138def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
139    """Setup a oblique grid of triangles
140    with m segments in the x-direction and n segments in the y-direction
141
142    """
143
144    from Numeric import array
145    import math
146
147    from config import epsilon
148
149
150    deltax = lenx/float(m)
151    deltay = leny/float(n)
152    a = 0.75*lenx*math.tan(theta/180.*math.pi)
153    x1 = lenx
154    y1 = 0
155    x2 = lenx
156    y2 = leny
157    x3 = 0.25*lenx
158    y3 = leny
159    x4 = x3
160    y4 = 0
161    a2 = a/(x1-x4)
162    a1 = -a2*x4
163    a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3)
164    a3 = 1. - (a1 + a2*x3)/y3 - a4*x3
165
166    # Dictionary of vertex objects
167    vertices = {}
168    points = []
169
170    for i in range(m+1):
171        x = deltax*i
172        for j in range(n+1):
173            y = deltay*j
174            if x > 0.25*lenx:
175                y = a1 + a2*x + a3*y + a4*x*y
176
177            vertices[i,j] = len(points)
178            points.append([x + origin[0], y + origin[1]])
179
180    # Construct 2 triangles per element
181    elements = []
182    boundary = {}
183    for i in range(m):
184        for j in range(n):
185            v1 = vertices[i,j+1]
186            v2 = vertices[i,j]
187            v3 = vertices[i+1,j+1]
188            v4 = vertices[i+1,j]
189
190            #Update boundary dictionary and create elements
191            if i == m-1:
192                boundary[(len(elements), 2)] = 'right'
193            if j == 0:
194                boundary[(len(elements), 1)] = 'bottom'
195            elements.append([v4,v3,v2]) #Lower
196
197            if i == 0:
198                boundary[(len(elements), 2)] = 'left'
199            if j == n-1:
200                boundary[(len(elements), 1)] = 'top'
201
202            elements.append([v1,v2,v3]) #Upper
203
204    return points, elements, boundary
205
206
207def contracting_channel(m, n, W_upstream = 1., W_downstream = 0.75,
208                            L_1 = 5.0, L_2 = 2.0, L_3 = 10, origin = (0.0, 0.0)):
209    """Setup a contracting channel grid of triangles
210    with m segments in the x-direction and n segments in the y-direction
211
212    """
213
214    from Numeric import array
215    import math
216
217    from config import epsilon
218
219
220    lenx = L_1 + L_2 + L_3
221    leny = W_upstream
222    deltax = lenx/float(m)
223    deltay = leny/float(n)
224   
225    x1 = 0
226    y1 = 0
227    x2 = L_1
228    y2 = 0
229    x3 = L_1 + L_2
230    y3 = (W_upstream - W_downstream)/2
231    x4 = L_1 + L_2 + L_3
232    y4 = y3
233    x5 = x4
234    y5 = y4 + W_downstream
235    x6 = L_1 + L_2
236    y6 = y5
237    x7 = L_1
238    y7 = W_upstream
239    x8 = 0
240    y8 = W_upstream
241    a1 = 0
242    a2 = (W_upstream - W_downstream)/(2*L_2)
243    a3 = 1
244    a4 = (W_downstream - W_upstream)/(L_2*W_upstream)
245   
246    # Dictionary of vertex objects
247    vertices = {}
248    points = []
249
250    for i in range(m+1):
251        x = deltax*i
252        for j in range(n+1):
253            y = deltay*j
254            if x > L_1 and x <= (L_1 + L_2):             
255                y = a1 + a2*(x - L_1) + a3*y + a4*(x - L_1)*y
256            elif x > L_1 + L_2:
257                y = (W_upstream - W_downstream)/2 + deltay*j*W_downstream/W_upstream
258
259            vertices[i,j] = len(points)
260            points.append([x + origin[0], y + origin[1]])
261
262    # Construct 2 triangles per element
263    elements = []
264    boundary = {}
265    for i in range(m):
266        for j in range(n):
267            v1 = vertices[i,j+1]
268            v2 = vertices[i,j]
269            v3 = vertices[i+1,j+1]
270            v4 = vertices[i+1,j]
271
272            #Update boundary dictionary and create elements
273            if i == m-1:
274                boundary[(len(elements), 2)] = 'right'
275            if j == 0:
276                boundary[(len(elements), 1)] = 'bottom'
277            elements.append([v4,v3,v2]) #Lower
278
279            if i == 0:
280                boundary[(len(elements), 2)] = 'left'
281            if j == n-1:
282                boundary[(len(elements), 1)] = 'top'
283
284            elements.append([v1,v2,v3]) #Upper
285
286    return points, elements, boundary
287
288
289
290def circular(m, n, radius=1.0, center = (0.0, 0.0)):
291    """Setup a circular grid of triangles with m concentric circles and
292    with n radial segments. If radius is are omitted the mesh defaults to
293    the unit circle radius.
294
295    radius: radius of circle
296
297    #FIXME: The triangles become degenerate for large values of m or n.
298    """
299
300
301
302    from math import pi, cos, sin
303
304    radius = float(radius) #Ensure floating point format
305
306    #Dictionary of vertex objects and list of points
307    vertices = {}
308    points = [[0.0, 0.0]] #Center point
309    vertices[0, 0] = 0
310
311    for i in range(n):
312        theta = 2*i*pi/n
313        x = cos(theta)
314        y = sin(theta)
315        for j in range(1,m+1):
316            delta = j*radius/m
317            vertices[i,j] = len(points)
318            points.append([delta*x, delta*y])
319
320    #Construct 2 triangles per element
321    elements = []
322    for i in range(n):
323        for j in range(1,m):
324
325            i1 = (i + 1) % n  #Wrap around
326
327            v1 = vertices[i,j+1]
328            v2 = vertices[i,j]
329            v3 = vertices[i1,j+1]
330            v4 = vertices[i1,j]
331
332            elements.append([v4,v2,v3]) #Lower
333            elements.append([v1,v3,v2]) #Upper
334
335
336    #Do the center
337    v1 = vertices[0,0]
338    for i in range(n):
339        i1 = (i + 1) % n  #Wrap around
340        v2 = vertices[i,1]
341        v3 = vertices[i1,1]
342
343        elements.append([v1,v2,v3]) #center
344
345    return points, elements
346
347
348def from_polyfile(name):
349    """Read mesh from .poly file, an obj like file format
350    listing first vertex coordinates and then connectivity
351    """
352
353    from util import anglediff
354    from math import pi
355    import os.path
356    root, ext = os.path.splitext(name)
357
358    if ext == 'poly':
359        filename = name
360    else:
361        filename = name + '.poly'
362
363
364    fid = open(filename)
365
366    points = []    #x, y
367    values = []    #z
368    ##vertex_values = []    #Repeated z
369    triangles = [] #v0, v1, v2
370
371    lines = fid.readlines()
372
373    keyword = lines[0].strip()
374    msg = 'First line in .poly file must contain the keyword: POINTS'
375    assert keyword == 'POINTS', msg
376
377    offending = 0
378    i = 1
379    while keyword == 'POINTS':
380        line = lines[i].strip()
381        i += 1
382
383        if line == 'POLYS':
384            keyword = line
385            break
386
387        fields = line.split(':')
388        assert int(fields[0]) == i-1, 'Point indices not consecutive'
389
390        #Split the three floats
391        xyz = fields[1].split()
392
393        x = float(xyz[0])
394        y = float(xyz[1])
395        z = float(xyz[2])
396
397        points.append([x, y])
398        values.append(z)
399
400
401    k = i
402    while keyword == 'POLYS':
403        line = lines[i].strip()
404        i += 1
405
406        if line == 'END':
407            keyword = line
408            break
409
410
411        fields = line.split(':')
412        assert int(fields[0]) == i-k, 'Poly indices not consecutive'
413
414        #Split the three indices
415        vvv = fields[1].split()
416
417        i0 = int(vvv[0])-1
418        i1 = int(vvv[1])-1
419        i2 = int(vvv[2])-1
420
421        #Check for and exclude degenerate areas
422        x0 = points[i0][0]
423        y0 = points[i0][1]
424        x1 = points[i1][0]
425        y1 = points[i1][1]
426        x2 = points[i2][0]
427        y2 = points[i2][1]
428
429        area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
430        if area > 0:
431
432            #Ensure that points are arranged in counter clock-wise order
433            v0 = [x1-x0, y1-y0]
434            v1 = [x2-x1, y2-y1]
435            v2 = [x0-x2, y0-y2]
436
437            a0 = anglediff(v1, v0)
438            a1 = anglediff(v2, v1)
439            a2 = anglediff(v0, v2)
440
441
442            if a0 < pi and a1 < pi and a2 < pi:
443                #all is well
444                j0 = i0
445                j1 = i1
446                j2 = i2
447            else:
448                #Swap two vertices
449                j0 = i1
450                j1 = i0
451                j2 = i2
452
453            triangles.append([j0, j1, j2])
454            ##vertex_values.append([values[j0], values[j1], values[j2]])
455        else:
456            offending +=1
457
458    print 'Removed %d offending triangles out of %d' %(offending, len(lines))
459    return points, triangles, values
460
461
462
463def strang_mesh(filename):
464    """Read Strang generated mesh.
465    """
466
467    from math import pi
468    from util import anglediff
469
470
471    fid = open(filename)
472    points = []    # List of x, y coordinates
473    triangles = [] # List of vertex ids as listed in the file
474
475    for line in fid.readlines():
476        fields = line.split()
477        if len(fields) == 2:
478            # we are reading vertex coordinates
479            points.append([float(fields[0]), float(fields[1])])
480        elif len(fields) == 3:
481            # we are reading triangle point id's (format ae+b)
482            triangles.append([int(float(fields[0]))-1,
483                              int(float(fields[1]))-1,
484                              int(float(fields[2]))-1])
485        else:
486            raise 'wrong format in ' + filename
487
488    elements = [] #Final list of elements
489
490    for t in triangles:
491        #Get vertex coordinates
492        v0 = t[0]
493        v1 = t[1]
494        v2 = t[2]
495
496        x0 = points[v0][0]
497        y0 = points[v0][1]
498        x1 = points[v1][0]
499        y1 = points[v1][1]
500        x2 = points[v2][0]
501        y2 = points[v2][1]
502
503        #Check that points are arranged in counter clock-wise order
504        vec0 = [x1-x0, y1-y0]
505        vec1 = [x2-x1, y2-y1]
506        vec2 = [x0-x2, y0-y2]
507
508        a0 = anglediff(vec1, vec0)
509        a1 = anglediff(vec2, vec1)
510        a2 = anglediff(vec0, vec2)
511
512        if a0 < pi and a1 < pi and a2 < pi:
513            elements.append([v0, v1, v2])
514        else:
515            elements.append([v0, v2, v1])
516
517    return points, elements
518
519# #Map from edge number to indices of associated vertices
520# edge_map = ((1,2), (0,2), (0,1))
521
522
Note: See TracBrowser for help on using the repository browser.