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

Last change on this file since 4901 was 4901, checked in by duncan, 16 years ago

Using a numeric array to store mesh triangulation info. The old way is also in place.

File size: 7.2 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
13from Numeric import array, Float, Int32, reshape
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, 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    try:
61        # If Int is used, instead of Int32, it fails in Linux
62        segments = ensure_numeric(segments, Int32)
63       
64    except ValueError:
65        msg = 'ERROR: Inconsistent segments array.'
66        raise ANUGAError, msg
67   
68    # This is after segments is numeric
69    if segatts is None or segatts == []:
70        segatts = [0 for x in range(segments.shape[0])]
71    try:
72        holes = ensure_numeric(holes, Float)
73    except ValueError:
74        msg = 'ERROR: Inconsistent holess array.'
75        raise ANUGAError, msg
76
77   
78    regions = add_area_tag(regions)
79    try:
80        regions = ensure_numeric(regions, Float)
81    except  (ValueError, TypeError):
82        msg = 'ERROR: Inconsistent regions array.'
83        raise ANUGAError, msg
84       
85    if not regions.shape[0] == 0 and regions.shape[1] <= 2:
86        msg = 'ERROR: Bad shape points array.'
87        raise ANUGAError, msg
88   
89    try:
90        pointatts = ensure_numeric(pointatts, Float)
91    except (ValueError, TypeError):
92        msg = 'ERROR: Inconsistent point attributes array.'
93        raise ANUGAError, msg
94
95    if pointatts.shape[0] <> points.shape[0]:
96        msg = """ERROR: Point attributes array not the same shape as
97        point array."""
98        raise ANUGAError, msg
99    if len(pointatts.shape) == 1:
100        pointatts = reshape(pointatts,(pointatts.shape[0],1))
101   
102    try:
103        segatts = ensure_numeric(segatts, Int32)
104    except ValueError:
105        msg = 'ERROR: Inconsistent point attributes array.'
106        raise ANUGAError, msg
107    if segatts.shape[0] <> segments.shape[0]:
108        msg = """ERROR: Segment attributes array not the same shape as
109        segment array."""
110        raise ANUGAError, msg
111
112   
113    #print "mode", mode
114    if mode.find('n'):
115        #pass
116        mode = 'j' + mode
117        # j- Jettisons vertices that are not part of the final
118        #    triangulation from the output .node file (including duplicate
119        #    input vertices and vertices ``eaten'' by holes).  - output a
120        #    list of neighboring triangles
121        # EG handles lone verts!
122           
123    #print "points",points
124    #print "segments", segments
125    #print "segments.shape", segments.shape
126    #print "holes", holes
127    #print "regions", regions
128    #print "pointatts", pointatts
129    #print "segatts", segatts
130    #XSprint "mode", mode
131    #print "yeah"
132    trianglelist, pointlist, pointmarkerlist, pointattributelist, triangleattributelist, segmentlist, segmentmarkerlist, neighborlist = triang.genMesh(points,segments,holes,regions,
133                          pointatts,segatts, mode, segments.flat)
134    mesh_dict = {}
135    # the values as arrays
136    mesh_dict['generatedtrianglelist'] = trianglelist
137    mesh_dict['generatedpointlist'] = pointlist
138    #print "mesh engine mesh_dict['generatedpointlist']", mesh_dict['generatedpointlist']
139    # WARNING - generatedpointmarkerlist IS UNTESTED
140    mesh_dict['generatedpointmarkerlist'] = pointmarkerlist
141    mesh_dict['generatedpointattributelist'] = pointattributelist
142    mesh_dict['generatedsegmentlist'] = segmentlist
143    mesh_dict['generatedsegmentmarkerlist'] =  segmentmarkerlist
144    mesh_dict['generatedtriangleneighborlist'] = neighborlist
145    mesh_dict['qaz'] = 1 #debugging
146
147    mesh_dict['triangleattributelist'] = triangleattributelist
148    #print "mesh eng generatedtrianglelist", trianglelist
149    #print "mesh eng mesh_dict['triangleattributelist'] ",mesh_dict['triangleattributelist']
150    #print "mesh eng mesh_dict['generatedtriangleattributelist'] ", mesh_dict['generatedtriangleattributelist']   
151    if True: 
152        mesh_dict['generatedtriangleattributelist'] = triangleattributelist
153
154        if mesh_dict['generatedtriangleattributelist'].shape[1] == 0:
155            mesh_dict['generatedtriangleattributelist'] = None
156           
157        if trianglelist.shape[0] == 0:
158            # There are no triangles.
159            # this is used by urs_ungridded2sww
160            raise NoTrianglesError
161    #print "mesh eng mesh_dict['generatedtriangleattributelist'] ", mesh_dict['generatedtriangleattributelist'] 
162   
163    # the structure of generatedtriangleattributelist is an list of
164    # list of integers.  It is transformed into a list of list of
165    # strings later on.  This is then inputted into an triangle
166    # object.  The triangle object outputs a list of strings.  Note
167    # the subtle change!  How should I handle this?  For my code, when
168    # the list of list of integers is transformed, transform it into a
169    # list of strings, not a list of list of strings.
170   
171    return mesh_dict
172
173def add_area_tag(regions):
174    """
175    So, what is the format?
176    A list with
177    [x,y,region_tag,area] OR [x,y,region_tag]
178    if it's [x,y,region_tag], add a 4th element, value of 0.0.
179    """
180    if isinstance(regions, ListType):
181        for i, region in enumerate(regions):
182            if len(region) == 3:
183                if isinstance(region, TupleType):
184                    #FIXME: How do you convert a tuple to a list?
185                    # I can do it a stupid way..
186                    tuple = region[:]
187                    regions[i] = []
188                    for j in tuple:
189                        regions[i].append(j)
190                    regions[i].append(0.0)
191                else:
192                    regions[i].append(0.0)
193                   
194                # let ensure numeric catch this..   
195                #len(region) <= 2:
196                #msg = 'ERROR: Inconsistent regions array.'
197                #raise msg
198            #elif
199    return regions
200
201if __name__ == "__main__":
202    pass 
Note: See TracBrowser for help on using the repository browser.