source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py @ 8124

Last change on this file since 8124 was 8124, checked in by wilsonr, 13 years ago

Made changes described in ticket 359.

File size: 24.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
206
207def rectangular_cross_slit(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
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.
213
214    len1: x direction (left to right)
215    len2: y direction (bottom to top)
216
217    Return to lists: points and elements suitable for creating a Mesh or
218    FVMesh object, e.g. Mesh(points, elements)
219    """
220
221    from anuga.config import epsilon
222
223    delta1 = float(len1)/m
224    delta2 = float(len2)/n
225
226    #Dictionary of vertex objects
227    vertices = {}
228    points = []
229
230    for i in range(m+1):
231        for j in range(n+1):
232            vertices[i,j] = len(points)
233            points.append([delta1*i + origin[0], delta2*+ origin[1]])
234
235    # Construct 4 triangles per element
236    elements = []
237    boundary = {}
238    for i in range(m):
239        for j in range(n):
240            v1 = vertices[i,j+1]
241            v2 = vertices[i,j]
242            v3 = vertices[i+1,j+1]
243            v4 = vertices[i+1,j]
244            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
245            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
246
247            # Create centre point
248            v5 = len(points)
249            points.append([x, y])
250
251            #Create left triangle
252            if i == 0:
253                boundary[(len(elements), 1)] = 'left'
254            elements.append([v2,v5,v1])
255
256            #Create bottom triangle
257            if j == 0:
258                boundary[(len(elements), 1)] = 'bottom'
259            elements.append([v4,v5,v2])
260
261            #Create right triangle
262            if i == m-1:
263                boundary[(len(elements), 1)] = 'right'
264            elements.append([v3,v5,v4])
265
266            #Create top triangle
267            if j == n-1:
268                boundary[(len(elements), 1)] = 'top'
269            elements.append([v1,v5,v3])
270
271
272    return points, elements, boundary
273
274
275
276def rectangular_periodic(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)):
277
278
279    """Setup a rectangular grid of triangles
280    with m+1 by n+1 grid points
281    and side lengths len1, len2. If side lengths are omitted
282    the mesh defaults to the unit square, divided between all the
283    processors
284
285    len1: x direction (left to right)
286    len2: y direction (bottom to top)
287
288    """
289
290    processor = 0
291    numproc   = 1
292
293   
294    n = n_g
295    m_low  = -1
296    m_high = m_g +1
297
298    m = m_high - m_low
299
300    delta1 = float(len1_g)/m_g
301    delta2 = float(len2_g)/n_g
302
303    len1 = len1_g*float(m)/float(m_g)
304    len2 = len2_g
305    origin = ( origin_g[0]+float(m_low)/float(m_g)*len1_g, origin_g[1] )
306
307    #Calculate number of points
308    Np = (m+1)*(n+1)
309
310    class VIndex:
311
312        def __init__(self, n,m):
313            self.n = n
314            self.m = m
315
316        def __call__(self, i,j):
317            return j+i*(self.n+1)
318
319    class EIndex:
320
321        def __init__(self, n,m):
322            self.n = n
323            self.m = m
324
325        def __call__(self, i,j):
326            return 2*(j+i*self.n)
327
328
329    I = VIndex(n,m)
330    E = EIndex(n,m)
331
332    points = num.zeros( (Np,2), num.float)
333
334    for i in range(m+1):
335        for j in range(n+1):
336
337            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
338
339    #Construct 2 triangles per rectangular element and assign tags to boundary
340    #Calculate number of triangles
341    Nt = 2*m*n
342
343
344    elements = num.zeros( (Nt,3), num.int)
345    boundary = {}
346    Idgl = []
347    Idfl = []
348    Idgr = []
349    Idfr = []
350
351    full_send_dict = {}
352    ghost_recv_dict = {}
353    nt = -1
354    for i in range(m):
355        for j in range(n):
356
357            i1 = I(i,j+1)
358            i2 = I(i,j)
359            i3 = I(i+1,j+1)
360            i4 = I(i+1,j)
361
362            #Lower Element
363            nt = E(i,j)
364            if i == 0:
365                Idgl.append(nt)
366
367            if i == 1:
368                Idfl.append(nt)
369
370            if i == m-2:
371                Idfr.append(nt)
372
373            if i == m-1:
374                Idgr.append(nt)
375
376            if i == m-1:
377                if processor == numproc-1:
378                    boundary[nt, 2] = 'right'
379                else:
380                    boundary[nt, 2] = 'ghost'
381       
382            if j == 0:
383                boundary[nt, 1] = 'bottom'
384            elements[nt,:] = [i4,i3,i2]
385
386            #Upper Element
387            nt = E(i,j)+1
388            if i == 0:
389                Idgl.append(nt)
390
391            if i == 1:
392                Idfl.append(nt)
393
394            if i == m-2:
395                Idfr.append(nt)
396
397            if i == m-1:
398                Idgr.append(nt)
399
400            if i == 0:
401                if processor == 0:
402                    boundary[nt, 2] = 'left'
403                else:
404                    boundary[nt, 2] = 'ghost'
405            if j == n-1:
406                boundary[nt, 1] = 'top'
407            elements[nt,:] = [i1,i2,i3]
408
409    Idfl.extend(Idfr)
410    Idgr.extend(Idgl)
411
412    Idfl = num.array(Idfl, num.int)
413    Idgr = num.array(Idgr, num.int)
414   
415    full_send_dict[processor]  = [Idfl, Idfl]
416    ghost_recv_dict[processor] = [Idgr, Idgr]
417
418
419    return  points, elements, boundary, full_send_dict, ghost_recv_dict
420
421
422def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
423    """Setup a oblique grid of triangles
424    with m segments in the x-direction and n segments in the y-direction
425
426    """
427
428    import math
429
430    from anuga.config import epsilon
431
432
433    deltax = lenx/float(m)
434    deltay = leny/float(n)
435    a = 0.75*lenx*math.tan(theta/180.*math.pi)
436    x1 = lenx
437    y1 = 0
438    x2 = lenx
439    y2 = leny
440    x3 = 0.25*lenx
441    y3 = leny
442    x4 = x3
443    y4 = 0
444    a2 = a/(x1-x4)
445    a1 = -a2*x4
446    a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3)
447    a3 = 1. - (a1 + a2*x3)/y3 - a4*x3
448
449    # Dictionary of vertex objects
450    vertices = {}
451    points = []
452
453    for i in range(m+1):
454        x = deltax*i
455        for j in range(n+1):
456            y = deltay*j
457            if x > 0.25*lenx:
458                y = a1 + a2*x + a3*y + a4*x*y
459
460            vertices[i,j] = len(points)
461            points.append([x + origin[0], y + origin[1]])
462
463    # Construct 2 triangles per element
464    elements = []
465    boundary = {}
466    for i in range(m):
467        for j in range(n):
468            v1 = vertices[i,j+1]
469            v2 = vertices[i,j]
470            v3 = vertices[i+1,j+1]
471            v4 = vertices[i+1,j]
472
473            #Update boundary dictionary and create elements
474            if i == m-1:
475                boundary[(len(elements), 2)] = 'right'
476            if j == 0:
477                boundary[(len(elements), 1)] = 'bottom'
478            elements.append([v4,v3,v2]) #Lower
479
480            if i == 0:
481                boundary[(len(elements), 2)] = 'left'
482            if j == n-1:
483                boundary[(len(elements), 1)] = 'top'
484
485            elements.append([v1,v2,v3]) #Upper
486
487    return points, elements, boundary
488
489
490def circular(m, n, radius=1.0, center = (0.0, 0.0)):
491    """Setup a circular grid of triangles with m concentric circles and
492    with n radial segments. If radius is are omitted the mesh defaults to
493    the unit circle radius.
494
495    radius: radius of circle
496
497    #FIXME: The triangles become degenerate for large values of m or n.
498    """
499
500
501
502    from math import pi, cos, sin
503
504    radius = float(radius) #Ensure floating point format
505
506    #Dictionary of vertex objects and list of points
507    vertices = {}
508    points = [[0.0, 0.0]] #Center point
509    vertices[0, 0] = 0
510
511    for i in range(n):
512        theta = 2*i*pi/n
513        x = cos(theta)
514        y = sin(theta)
515        for j in range(1,m+1):
516            delta = j*radius/m
517            vertices[i,j] = len(points)
518            points.append([delta*x, delta*y])
519
520    #Construct 2 triangles per element
521    elements = []
522    for i in range(n):
523        for j in range(1,m):
524
525            i1 = (i + 1) % n  #Wrap around
526
527            v1 = vertices[i,j+1]
528            v2 = vertices[i,j]
529            v3 = vertices[i1,j+1]
530            v4 = vertices[i1,j]
531
532            elements.append([v4,v2,v3]) #Lower
533            elements.append([v1,v3,v2]) #Upper
534
535
536    #Do the center
537    v1 = vertices[0,0]
538    for i in range(n):
539        i1 = (i + 1) % n  #Wrap around
540        v2 = vertices[i,1]
541        v3 = vertices[i1,1]
542
543        elements.append([v1,v2,v3]) #center
544
545    return points, elements
546
547
548def from_polyfile(name):
549    """Read mesh from .poly file, an obj like file format
550    listing first vertex coordinates and then connectivity
551    """
552
553    from anuga.utilities.numerical_tools import anglediff
554    from math import pi
555    import os.path
556    root, ext = os.path.splitext(name)
557
558    if ext == 'poly':
559        filename = name
560    else:
561        filename = name + '.poly'
562
563
564    fid = open(filename)
565
566    points = []    #x, y
567    values = []    #z
568    ##vertex_values = []    #Repeated z
569    triangles = [] #v0, v1, v2
570
571    lines = fid.readlines()
572
573    keyword = lines[0].strip()
574    msg = 'First line in .poly file must contain the keyword: POINTS'
575    assert keyword == 'POINTS', msg
576
577    offending = 0
578    i = 1
579    while keyword == 'POINTS':
580        line = lines[i].strip()
581        i += 1
582
583        if line == 'POLYS':
584            keyword = line
585            break
586
587        fields = line.split(':')
588        assert int(fields[0]) == i-1, 'Point indices not consecutive'
589
590        #Split the three floats
591        xyz = fields[1].split()
592
593        x = float(xyz[0])
594        y = float(xyz[1])
595        z = float(xyz[2])
596
597        points.append([x, y])
598        values.append(z)
599
600
601    k = i
602    while keyword == 'POLYS':
603        line = lines[i].strip()
604        i += 1
605
606        if line == 'END':
607            keyword = line
608            break
609
610
611        fields = line.split(':')
612        assert int(fields[0]) == i-k, 'Poly indices not consecutive'
613
614        #Split the three indices
615        vvv = fields[1].split()
616
617        i0 = int(vvv[0])-1
618        i1 = int(vvv[1])-1
619        i2 = int(vvv[2])-1
620
621        #Check for and exclude degenerate areas
622        x0 = points[i0][0]
623        y0 = points[i0][1]
624        x1 = points[i1][0]
625        y1 = points[i1][1]
626        x2 = points[i2][0]
627        y2 = points[i2][1]
628
629        area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
630        if area > 0:
631
632            #Ensure that points are arranged in counter clock-wise order
633            v0 = [x1-x0, y1-y0]
634            v1 = [x2-x1, y2-y1]
635            v2 = [x0-x2, y0-y2]
636
637            a0 = anglediff(v1, v0)
638            a1 = anglediff(v2, v1)
639            a2 = anglediff(v0, v2)
640
641
642            if a0 < pi and a1 < pi and a2 < pi:
643                #all is well
644                j0 = i0
645                j1 = i1
646                j2 = i2
647            else:
648                #Swap two vertices
649                j0 = i1
650                j1 = i0
651                j2 = i2
652
653            triangles.append([j0, j1, j2])
654            ##vertex_values.append([values[j0], values[j1], values[j2]])
655        else:
656            offending +=1
657
658    log.critical('Removed %d offending triangles out of %d'
659                 % (offending, len(lines)))
660    return points, triangles, values
661
662
663
664def strang_mesh(filename):
665    """Read Strang generated mesh.
666    """
667
668    from math import pi
669    from anuga.utilities.numerical_tools import anglediff
670
671
672    fid = open(filename)
673    points = []    # List of x, y coordinates
674    triangles = [] # List of vertex ids as listed in the file
675
676    for line in fid.readlines():
677        fields = line.split()
678        if len(fields) == 2:
679            # we are reading vertex coordinates
680            points.append([float(fields[0]), float(fields[1])])
681        elif len(fields) == 3:
682            # we are reading triangle point id's (format ae+b)
683            triangles.append([int(float(fields[0]))-1,
684                              int(float(fields[1]))-1,
685                              int(float(fields[2]))-1])
686        else:
687            raise Excetion('wrong format in %s' % filename)
688
689    elements = [] #Final list of elements
690
691    for t in triangles:
692        #Get vertex coordinates
693        v0 = t[0]
694        v1 = t[1]
695        v2 = t[2]
696
697        x0 = points[v0][0]
698        y0 = points[v0][1]
699        x1 = points[v1][0]
700        y1 = points[v1][1]
701        x2 = points[v2][0]
702        y2 = points[v2][1]
703
704        #Check that points are arranged in counter clock-wise order
705        vec0 = [x1-x0, y1-y0]
706        vec1 = [x2-x1, y2-y1]
707        vec2 = [x0-x2, y0-y2]
708
709        a0 = anglediff(vec1, vec0)
710        a1 = anglediff(vec2, vec1)
711        a2 = anglediff(vec0, vec2)
712
713        if a0 < pi and a1 < pi and a2 < pi:
714            elements.append([v0, v1, v2])
715        else:
716            elements.append([v0, v2, v1])
717
718    return points, elements
719
720# #Map from edge number to indices of associated vertices
721# edge_map = ((1,2), (0,2), (0,1))
722
723def contracting_channel(m, n, W_upstream = 1., W_downstream = 0.75,
724                            L_1 = 5.0, L_2 = 2.0, L_3 = 10, origin = (0.0, 0.0)):
725    """Setup a contracting channel grid of triangles
726    with m segments in the x-direction and n segments in the y-direction
727
728    """
729
730    import math
731
732    from anuga.config import epsilon
733
734
735    lenx = L_1 + L_2 + L_3
736    leny = W_upstream
737    deltax = lenx/float(m)
738    deltay = leny/float(n)
739
740    x1 = 0
741    y1 = 0
742    x2 = L_1
743    y2 = 0
744    x3 = L_1 + L_2
745    y3 = (W_upstream - W_downstream)/2
746    x4 = L_1 + L_2 + L_3
747    y4 = y3
748    x5 = x4
749    y5 = y4 + W_downstream
750    x6 = L_1 + L_2
751    y6 = y5
752    x7 = L_1
753    y7 = W_upstream
754    x8 = 0
755    y8 = W_upstream
756    a1 = 0
757    a2 = (W_upstream - W_downstream)/(2*L_2)
758    a3 = 1
759    a4 = (W_downstream - W_upstream)/(L_2*W_upstream)
760
761    # Dictionary of vertex objects
762    vertices = {}
763    points = []
764
765    for i in range(m+1):
766        x = deltax*i
767        for j in range(n+1):
768            y = deltay*j
769            if x > L_1 and x <= (L_1 + L_2):
770                y = a1 + a2*(x - L_1) + a3*y + a4*(x - L_1)*y
771            elif x > L_1 + L_2:
772                y = (W_upstream - W_downstream)/2 + deltay*j*W_downstream/W_upstream
773
774            vertices[i,j] = len(points)
775            points.append([x + origin[0], y + origin[1]])
776
777    # Construct 2 triangles per element
778    elements = []
779    boundary = {}
780    for i in range(m):
781        for j in range(n):
782            v1 = vertices[i,j+1]
783            v2 = vertices[i,j]
784            v3 = vertices[i+1,j+1]
785            v4 = vertices[i+1,j]
786
787            #Update boundary dictionary and create elements
788            if i == m-1:
789                boundary[(len(elements), 2)] = 'right'
790            if j == 0:
791                boundary[(len(elements), 1)] = 'bottom'
792            elements.append([v4,v3,v2]) #Lower
793
794            if i == 0:
795                boundary[(len(elements), 2)] = 'left'
796            if j == n-1:
797                boundary[(len(elements), 1)] = 'top'
798
799            elements.append([v1,v2,v3]) #Upper
800
801    return points, elements, boundary
802
803
804def contracting_channel_cross(m, n, W_upstream = 1., W_downstream = 0.75,
805                              L_1 = 5.0, L_2 = 2.0, L_3 = 10, origin = (0.0, 0.0)):
806    """Setup a contracting channel grid of triangles
807    with m segments in the x-direction and n segments in the y-direction
808
809    """
810
811    import math
812
813    from anuga.config import epsilon
814
815
816    lenx = L_1 + L_2 + L_3
817    leny = W_upstream
818    deltax = lenx/float(m)
819    deltay = leny/float(n)
820
821    x1 = 0
822    y1 = 0
823    x2 = L_1
824    y2 = 0
825    x3 = L_1 + L_2
826    y3 = (W_upstream - W_downstream)/2
827    x4 = L_1 + L_2 + L_3
828    y4 = y3
829    x5 = x4
830    y5 = y4 + W_downstream
831    x6 = L_1 + L_2
832    y6 = y5
833    x7 = L_1
834    y7 = W_upstream
835    x8 = 0
836    y8 = W_upstream
837    a1 = 0
838    a2 = (W_upstream - W_downstream)/(2*L_2)
839    a3 = 1
840    a4 = (W_downstream - W_upstream)/(L_2*W_upstream)
841
842    # Dictionary of vertex objects
843    vertices = {}
844    points = []
845
846    for i in range(m+1):
847        x = deltax*i
848        for j in range(n+1):
849            y = deltay*j
850            if x > L_1 and x <= (L_1 + L_2):
851                y = a1 + a2*(x - L_1) + a3*y + a4*(x - L_1)*y
852            elif x > L_1 + L_2:
853                y = (W_upstream - W_downstream)/2 + deltay*j*W_downstream/W_upstream
854
855            vertices[i,j] = len(points)
856            points.append([x + origin[0], y + origin[1]])
857
858    # Construct 4 triangles per element
859    elements = []
860    boundary = {}
861    for i in range(m):
862        for j in range(n):
863            v1 = vertices[i,j+1]
864            v2 = vertices[i,j]
865            v3 = vertices[i+1,j+1]
866            v4 = vertices[i+1,j]
867            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
868            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
869            v5 = len(points)
870            points.append([x, y])
871
872            #Create left triangle
873            if i == 0:
874                boundary[(len(elements), 1)] = 'left'
875            elements.append([v2,v5,v1])
876
877            #Create bottom triangle
878            if j == 0:
879                boundary[(len(elements), 1)] = 'bottom'
880            elements.append([v4,v5,v2])
881
882            #Create right triangle
883            if i == m-1:
884                boundary[(len(elements), 1)] = 'right'
885            elements.append([v3,v5,v4])
886
887            #Create top triangle
888            if j == n-1:
889                boundary[(len(elements), 1)] = 'top'
890            elements.append([v1,v5,v3])
891
892
893    return points, elements, boundary
894
895
896
897
898def oblique_cross(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
899    """Setup a oblique grid of triangles
900    with m segments in the x-direction and n segments in the y-direction
901
902    """
903
904    import math
905
906    from anuga.config import epsilon
907
908
909    deltax = lenx/float(m)
910    deltay = leny/float(n)
911    a = 0.75*lenx*math.tan(theta/180.*math.pi)
912    x1 = lenx
913    y1 = 0
914    x2 = lenx
915    y2 = leny
916    x3 = 0.25*lenx
917    y3 = leny
918    x4 = x3
919    y4 = 0
920    a2 = a/(x1-x4)
921    a1 = -a2*x4
922    a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3)
923    a3 = 1. - (a1 + a2*x3)/y3 - a4*x3
924
925    # Dictionary of vertex objects
926    vertices = {}
927    points = []
928
929    for i in range(m+1):
930        x = deltax*i
931        for j in range(n+1):
932            y = deltay*j
933            if x > 0.25*lenx:
934                y = a1 + a2*x + a3*y + a4*x*y
935
936            vertices[i,j] = len(points)
937            points.append([x + origin[0], y + origin[1]])
938
939    # Construct 4 triangles per element
940    elements = []
941    boundary = {}
942    for i in range(m):
943        for j in range(n):
944            v1 = vertices[i,j+1]
945            v2 = vertices[i,j]
946            v3 = vertices[i+1,j+1]
947            v4 = vertices[i+1,j]
948            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
949            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
950            v5 = len(points)
951            points.append([x, y])
952
953            #Update boundary dictionary and create elements
954                        #Create left triangle
955            if i == 0:
956                boundary[(len(elements), 1)] = 'left'
957            elements.append([v2,v5,v1])
958
959            #Create bottom triangle
960            if j == 0:
961                boundary[(len(elements), 1)] = 'bottom'
962            elements.append([v4,v5,v2])
963
964            #Create right triangle
965            if i == m-1:
966                boundary[(len(elements), 1)] = 'right'
967            elements.append([v3,v5,v4])
968
969            #Create top triangle
970            if j == n-1:
971                boundary[(len(elements), 1)] = 'top'
972            elements.append([v1,v5,v3])
973
974
975    return points, elements, boundary
Note: See TracBrowser for help on using the repository browser.