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

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

Replaced 'print' statements with log.critical() calls.

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