source: anuga_core/source/anuga/mesh_engine/mesh_engine.py @ 7276

Last change on this file since 7276 was 7276, checked in by ole, 15 years ago

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

The suites test_all.py (in source/anuga) and validate_all.py passed using Python2.5 with numpy on my Ubuntu Linux box.

File size: 7.4 KB
Line 
1#!/usr/bin/env python
2
3import sys
4
5from types import ListType, TupleType
6
7import exceptions
8
9class NoTrianglesError(exceptions.Exception): pass
10import anuga.mesh_engine.mesh_engine_c_layer as triang
11#import anuga.mesh_engine.list_dic as triang
12
13import numpy as num
14
15from anuga.utilities.numerical_tools import ensure_numeric
16from anuga.utilities.anuga_exceptions import ANUGAError
17   
18def generate_mesh(points=None,
19                  segments=None,holes=None,regions=None,
20                  pointatts=None,segatts=None,
21                  mode=None, dummy_test=None):
22    """
23    pointatts can be a list of lists.
24
25    generatedtriangleattributelist is used to represent tagged regions.
26    #FIXME (DSG-DSG): add comments
27    """
28    #FIXME (DSG-DSG): Catch parameters that are lists,
29    #instead of lists of lists
30    # check shape[1] is 2 etc
31
32    if points is None:
33        points = []
34
35    if segments is None:
36        segments = []
37
38    if holes is None:
39        holes = []
40       
41    if regions is None:
42        regions = []
43
44    if dummy_test is None:
45        dummy_test  = []
46       
47    try:
48        points =  ensure_numeric(points, num.float)
49    except ValueError:
50        msg = 'ERROR: Inconsistent points array.'
51        raise ANUGAError, msg
52    if points.shape[1] <>2:
53        msg = 'ERROR: Bad shape points array.'
54        raise ANUGAError, msg
55
56    #print "pointatts",pointatts
57    # This is after points is numeric
58    if pointatts is None or pointatts == []:
59        pointatts = [[] for x in range(points.shape[0])]
60       
61    try:
62        # If num.int is used, instead of num.int32, it fails in Linux
63        segments = ensure_numeric(segments, num.int32)
64       
65    except ValueError:
66        msg = 'ERROR: Inconsistent segments array.'
67        raise ANUGAError, msg
68   
69    # This is after segments is numeric
70    if segatts is None or segatts == []:
71        segatts = [0 for x in range(segments.shape[0])]
72       
73    try:
74        holes = ensure_numeric(holes, num.float)
75    except ValueError:
76        msg = 'ERROR: Inconsistent holess array.'
77        raise ANUGAError, msg
78
79   
80    regions = add_area_tag(regions)
81    try:
82        regions = ensure_numeric(regions, num.float)
83    except  (ValueError, TypeError):
84        msg = 'ERROR: Inconsistent regions array.'
85        raise ANUGAError, msg
86       
87    if not regions.shape[0] == 0 and regions.shape[1] <= 2:
88        msg = 'ERROR: Bad shape points array.'
89        raise ANUGAError, msg
90   
91    try:
92        pointatts = ensure_numeric(pointatts, num.float)
93    except (ValueError, TypeError):
94        msg = 'ERROR: Inconsistent point attributes array.'
95        raise ANUGAError, msg
96
97    if pointatts.shape[0] <> points.shape[0]:
98        msg = """ERROR: Point attributes array not the same shape as
99        point array."""
100        raise ANUGAError, msg
101    if len(pointatts.shape) == 1:
102        pointatts = num.reshape(pointatts,(pointatts.shape[0],1))
103   
104    try:
105        segatts = ensure_numeric(segatts, num.int32)
106    except ValueError:
107        msg = 'ERROR: Inconsistent point attributes array.'
108        raise ANUGAError, msg
109    if segatts.shape[0] <> segments.shape[0]:
110        msg = """ERROR: Segment attributes array not the same shape as
111        segment array."""
112        raise ANUGAError, msg
113   
114    #print "mode", mode
115    if mode.find('n'):
116        #pass
117        mode = 'j' + mode
118        # j- Jettisons vertices that are not part of the final
119        #    triangulation from the output .node file (including duplicate
120        #    input vertices and vertices ``eaten'' by holes).  - output a
121        #    list of neighboring triangles
122        # EG handles lone verts!
123           
124    #print "points",points
125    #print "segments", segments
126    #print "segments.shape", segments.shape
127    #print "holes", holes
128    #print "regions", regions
129    #print "pointatts", pointatts
130    #print "segatts", segatts
131    #print "mode", mode
132    #print "yeah"
133    # .copy()
134    trianglelist, pointlist, pointmarkerlist, pointattributelist, triangleattributelist, segmentlist, segmentmarkerlist, neighborlist = triang.genMesh(points,segments,holes,regions,
135                          pointatts,segatts, mode)
136    mesh_dict = {}
137    # the values as arrays
138    mesh_dict['generatedtrianglelist'] = trianglelist
139    mesh_dict['generatedpointlist'] = pointlist
140    #print "mesh engine mesh_dict['generatedpointlist']", mesh_dict['generatedpointlist']
141    # WARNING - generatedpointmarkerlist IS UNTESTED
142    mesh_dict['generatedpointmarkerlist'] = pointmarkerlist
143    mesh_dict['generatedpointattributelist'] = pointattributelist
144    mesh_dict['generatedsegmentlist'] = segmentlist
145    mesh_dict['generatedsegmentmarkerlist'] =  segmentmarkerlist
146    mesh_dict['generatedtriangleneighborlist'] = neighborlist
147    mesh_dict['qaz'] = 1 #debugging
148
149    #mesh_dict['triangleattributelist'] = triangleattributelist
150    #print "mesh eng generatedtrianglelist", trianglelist
151    #print "mesh eng mesh_dict['triangleattributelist'] ",mesh_dict['triangleattributelist']
152    #print "mesh eng mesh_dict['generatedtriangleattributelist'] ", mesh_dict['generatedtriangleattributelist']   
153    if True: 
154        mesh_dict['generatedtriangleattributelist'] = triangleattributelist
155
156        if mesh_dict['generatedtriangleattributelist'].shape[1] == 0:
157            mesh_dict['generatedtriangleattributelist'] = None
158           
159        if mesh_dict['generatedpointattributelist'].shape[1] == 0:
160            mesh_dict['generatedpointattributelist'] = None
161           
162        if mesh_dict['generatedtriangleneighborlist'].shape[1] == 0:
163            mesh_dict['generatedtriangleneighborlist'] = None
164           
165        if trianglelist.shape[0] == 0:
166            # There are no triangles.
167            # this is used by urs_ungridded2sww
168            raise NoTrianglesError
169    #print "mesh eng mesh_dict['generatedtriangleattributelist'] ", mesh_dict['generatedtriangleattributelist'] 
170    a = mesh_dict['generatedtriangleattributelist']
171    #print 'mesh_dict', mesh_dict
172    # the structure of generatedtriangleattributelist is an list of
173    # list of integers.  It is transformed into a list of list of
174    # strings later on.  This is then inputted into an triangle
175    # object.  The triangle object outputs a list of strings.  Note
176    # the subtle change!  How should I handle this?  For my code, when
177    # the list of list of integers is transformed, transform it into a
178    # list of strings, not a list of list of strings.
179   
180    return mesh_dict
181
182def add_area_tag(regions):
183    """
184    So, what is the format?
185    A list with
186    [x,y,region_tag,area] OR [x,y,region_tag]
187    if it's [x,y,region_tag], add a 4th element, value of 0.0.
188    """
189    if isinstance(regions, ListType):
190        for i, region in enumerate(regions):
191            if len(region) == 3:
192                if isinstance(region, TupleType):
193                    #FIXME: How do you convert a tuple to a list?
194                    # I can do it a stupid way..
195                    tuple = region[:]
196                    regions[i] = []
197                    for j in tuple:
198                        regions[i].append(j)
199                    regions[i].append(0.0)
200                else:
201                    regions[i].append(0.0)
202                   
203                # let ensure numeric catch this..   
204                #len(region) <= 2:
205                #msg = 'ERROR: Inconsistent regions array.'
206                #raise msg
207            #elif
208    return regions
209
210if __name__ == "__main__":
211    pass 
Note: See TracBrowser for help on using the repository browser.