source: inundation/ga/storm_surge/pyvolution/mesh_factory.py @ 1520

Last change on this file since 1520 was 1387, checked in by steve, 20 years ago

Need to look at uint test for inscribed circle

File size: 11.3 KB
RevLine 
[1387]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 circular(m, n, radius=1.0, center = (0.0, 0.0)):
208    """Setup a circular grid of triangles with m concentric circles and
209    with n radial segments. If radius is are omitted the mesh defaults to
210    the unit circle radius.
211
212    radius: radius of circle
213
214    #FIXME: The triangles become degenerate for large values of m or n.
215    """
216
217
218
219    from math import pi, cos, sin
220
221    radius = float(radius) #Ensure floating point format
222
223    #Dictionary of vertex objects and list of points
224    vertices = {}
225    points = [[0.0, 0.0]] #Center point
226    vertices[0, 0] = 0
227
228    for i in range(n):
229        theta = 2*i*pi/n
230        x = cos(theta)
231        y = sin(theta)
232        for j in range(1,m+1):
233            delta = j*radius/m
234            vertices[i,j] = len(points)
235            points.append([delta*x, delta*y])
236
237    #Construct 2 triangles per element
238    elements = []
239    for i in range(n):
240        for j in range(1,m):
241
242            i1 = (i + 1) % n  #Wrap around
243
244            v1 = vertices[i,j+1]
245            v2 = vertices[i,j]
246            v3 = vertices[i1,j+1]
247            v4 = vertices[i1,j]
248
249            elements.append([v4,v2,v3]) #Lower
250            elements.append([v1,v3,v2]) #Upper
251
252
253    #Do the center
254    v1 = vertices[0,0]
255    for i in range(n):
256        i1 = (i + 1) % n  #Wrap around
257        v2 = vertices[i,1]
258        v3 = vertices[i1,1]
259
260        elements.append([v1,v2,v3]) #center
261
262    return points, elements
263
264
265def from_polyfile(name):
266    """Read mesh from .poly file, an obj like file format
267    listing first vertex coordinates and then connectivity
268    """
269
270    from util import anglediff
271    from math import pi
272    import os.path
273    root, ext = os.path.splitext(name)
274
275    if ext == 'poly':
276        filename = name
277    else:
278        filename = name + '.poly'
279
280
281    fid = open(filename)
282
283    points = []    #x, y
284    values = []    #z
285    ##vertex_values = []    #Repeated z
286    triangles = [] #v0, v1, v2
287
288    lines = fid.readlines()
289
290    keyword = lines[0].strip()
291    msg = 'First line in .poly file must contain the keyword: POINTS'
292    assert keyword == 'POINTS', msg
293
294    offending = 0
295    i = 1
296    while keyword == 'POINTS':
297        line = lines[i].strip()
298        i += 1
299
300        if line == 'POLYS':
301            keyword = line
302            break
303
304        fields = line.split(':')
305        assert int(fields[0]) == i-1, 'Point indices not consecutive'
306
307        #Split the three floats
308        xyz = fields[1].split()
309
310        x = float(xyz[0])
311        y = float(xyz[1])
312        z = float(xyz[2])
313
314        points.append([x, y])
315        values.append(z)
316
317
318    k = i
319    while keyword == 'POLYS':
320        line = lines[i].strip()
321        i += 1
322
323        if line == 'END':
324            keyword = line
325            break
326
327
328        fields = line.split(':')
329        assert int(fields[0]) == i-k, 'Poly indices not consecutive'
330
331        #Split the three indices
332        vvv = fields[1].split()
333
334        i0 = int(vvv[0])-1
335        i1 = int(vvv[1])-1
336        i2 = int(vvv[2])-1
337
338        #Check for and exclude degenerate areas
339        x0 = points[i0][0]
340        y0 = points[i0][1]
341        x1 = points[i1][0]
342        y1 = points[i1][1]
343        x2 = points[i2][0]
344        y2 = points[i2][1]
345
346        area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
347        if area > 0:
348
349            #Ensure that points are arranged in counter clock-wise order
350            v0 = [x1-x0, y1-y0]
351            v1 = [x2-x1, y2-y1]
352            v2 = [x0-x2, y0-y2]
353
354            a0 = anglediff(v1, v0)
355            a1 = anglediff(v2, v1)
356            a2 = anglediff(v0, v2)
357
358
359            if a0 < pi and a1 < pi and a2 < pi:
360                #all is well
361                j0 = i0
362                j1 = i1
363                j2 = i2
364            else:
365                #Swap two vertices
366                j0 = i1
367                j1 = i0
368                j2 = i2
369
370            triangles.append([j0, j1, j2])
371            ##vertex_values.append([values[j0], values[j1], values[j2]])
372        else:
373            offending +=1
374
375    print 'Removed %d offending triangles out of %d' %(offending, len(lines))
376    return points, triangles, values
377
378
379
380def strang_mesh(filename):
381    """Read Strang generated mesh.
382    """
383
384    from math import pi
385    from util import anglediff
386
387
388    fid = open(filename)
389    points = []    # List of x, y coordinates
390    triangles = [] # List of vertex ids as listed in the file
391
392    for line in fid.readlines():
393        fields = line.split()
394        if len(fields) == 2:
395            # we are reading vertex coordinates
396            points.append([float(fields[0]), float(fields[1])])
397        elif len(fields) == 3:
398            # we are reading triangle point id's (format ae+b)
399            triangles.append([int(float(fields[0]))-1,
400                              int(float(fields[1]))-1,
401                              int(float(fields[2]))-1])
402        else:
403            raise 'wrong format in ' + filename
404
405    elements = [] #Final list of elements
406
407    for t in triangles:
408        #Get vertex coordinates
409        v0 = t[0]
410        v1 = t[1]
411        v2 = t[2]
412
413        x0 = points[v0][0]
414        y0 = points[v0][1]
415        x1 = points[v1][0]
416        y1 = points[v1][1]
417        x2 = points[v2][0]
418        y2 = points[v2][1]
419
420        #Check that points are arranged in counter clock-wise order
421        vec0 = [x1-x0, y1-y0]
422        vec1 = [x2-x1, y2-y1]
423        vec2 = [x0-x2, y0-y2]
424
425        a0 = anglediff(vec1, vec0)
426        a1 = anglediff(vec2, vec1)
427        a2 = anglediff(vec0, vec2)
428
429        if a0 < pi and a1 < pi and a2 < pi:
430            elements.append([v0, v1, v2])
431        else:
432            elements.append([v0, v2, v1])
433
434    return points, elements
435
436# #Map from edge number to indices of associated vertices
437# edge_map = ((1,2), (0,2), (0,1))
438
439
Note: See TracBrowser for help on using the repository browser.