source: inundation/pyvolution/pmesh2domain.py @ 2267

Last change on this file since 2267 was 2056, checked in by duncan, 19 years ago

create a domain instance from a pmesh.mesh instance

File size: 6.6 KB
RevLine 
[374]1"""Class pmesh2domain - Converting .tsh files to doamains
[371]2
3
4   Copyright 2004
5   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
6   Geoscience Australia
7"""
8
[2056]9# FXME (DSG-DSG): Combine with function pmesh_to_domain_instance
10def pmesh_instance_to_domain_instance(mesh,
11                                      DomainClass,
12                                      setting_function=None):
13    """
14    """
15    import sys
16    from domain import Domain
[371]17
[2056]18    vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \
19                        ,tagged_elements_dict, geo_reference = \
20                pmesh_to_domain(mesh_instance=mesh,
21                                setting_function=setting_function)
22
23
24    assert issubclass(DomainClass, Domain),"DomainClass is not a subclass of Domain."
25
26
27    domain = DomainClass(coordinates = vertex_coordinates,
28                         vertices = vertices,
29                         boundary = tag_dict,
30                         tagged_elements = tagged_elements_dict,
31                         geo_reference = geo_reference )
32
33    # set the water stage to be the elevation
34    if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'):
35        vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']
36
37    domain.set_quantity_vertices_dict(vertex_quantity_dict)
38    #print "vertex_quantity_dict",vertex_quantity_dict
39    return domain
40
[415]41def pmesh_to_domain_instance(fileName, DomainClass, setting_function = None):
[371]42    """
43    """
44    import sys
[389]45    from domain import Domain
[371]46
[1556]47    vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \
[1178]48                        ,tagged_elements_dict, geo_reference = \
[2056]49                pmesh_to_domain(fileName=fileName,
50                                setting_function=setting_function)
[1556]51
52
[1178]53    assert issubclass(DomainClass, Domain),"DomainClass is not a subclass of Domain."
[636]54
55
[1556]56    domain = DomainClass(coordinates = vertex_coordinates,
57                         vertices = vertices,
58                         boundary = tag_dict,
[1178]59                         tagged_elements = tagged_elements_dict,
60                         geo_reference = geo_reference )
[371]61
[636]62
[1556]63
[636]64    #FIXME (Ole): Is this really the right place to apply the a default
65    #value specific to the shallow water wave equation?
[1556]66    #The 'assert' above indicates that any subclass of Domain is acceptable.
[636]67    #Suggestion - module shallow_water.py will eventually take care of this
[1556]68    #(when I get around to it) so it should be removed from here.
69
[905]70    # This doesn't work on the domain instance.
71    # This is still needed so -ve elevations don't cuase 'lakes'
[1178]72    # The fixme we discussed was to only create a quantity when its values
73    #are set.
[905]74    # I think that's the way to go still
[1556]75
[773]76    # set the water stage to be the elevation
77    if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'):
78        vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']
[1556]79
[389]80    domain.set_quantity_vertices_dict(vertex_quantity_dict)
81    #print "vertex_quantity_dict",vertex_quantity_dict
82    return domain
[371]83
84
[389]85
[2056]86def pmesh_to_domain(fileName=None,
87                    mesh_instance=None,
88                    setting_function=None):
[371]89    """
90    convert a pmesh dictionary to a list of Volumes.
91    Also, return a list of triangles which have boundary tags
[389]92    mesh_dict structure;
[1556]93    generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
[371]94    generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
95    generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
[969]96    generated segment tag list: [S1Tag, S2Tag, ...] (list of ints)
[371]97    triangle list:  [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers)
98    triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor
[415]99    triangle attribute list: [T1att, T2att, ...] (list of strings)
[371]100    """
[1556]101
[389]102    from Numeric import transpose
[1423]103    from load_mesh.loadASCII import import_mesh_file
[1556]104
[2056]105    if fileName is None:
106        mesh_dict = mesh_instance.Mesh2IODict()
107    else:
108        mesh_dict = import_mesh_file(fileName)
[389]109    #print "mesh_dict",mesh_dict
[968]110    vertex_coordinates = mesh_dict['vertices']
111    volumes = mesh_dict['triangles']
[389]112    vertex_quantity_dict = {}
[968]113    point_atts = transpose(mesh_dict['vertex_attributes'])
114    point_titles  = mesh_dict['vertex_attribute_titles']
[1178]115    geo_reference  = mesh_dict['geo_reference']
[389]116    for quantity, value_vector in map (None, point_titles, point_atts):
117        vertex_quantity_dict[quantity] = value_vector
[969]118    tag_dict = pmesh_dict_to_tag_dict(mesh_dict)
[415]119    tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict)
[1556]120    return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, tagged_elements_dict, geo_reference
[371]121
122
123
[415]124def build_tagged_elements_dictionary(mesh_dict):
125    """Build the dictionary of element tags.
126    tagged_elements is a dictionary of element arrays,
127    keyed by tag:
128    { (tag): [e1, e2, e3..] }
129    """
[969]130    tri_atts = mesh_dict['triangle_tags']
[415]131    #print "tri_atts", tri_atts
132    tagged_elements = {}
133    for tri_att_index in range(len(tri_atts)):
[1183]134        tagged_elements.setdefault(tri_atts[tri_att_index],[]).append(tri_att_index)
[415]135    #print "DSG pm2do tagged_elements", tagged_elements
136    return tagged_elements
[371]137
[969]138def pmesh_dict_to_tag_dict(mesh_dict):
[389]139    """ Convert the pmesh dictionary (mesh_dict) description of boundary tags
[969]140    to a dictionary of tags, indexed with volume id and face number.
[371]141    """
[968]142    triangles = mesh_dict['triangles']
[371]143    sides = calc_sides(triangles)
[969]144    tag_dict = {}
145    for seg, tag in map(None,mesh_dict['segments'],
146                           mesh_dict['segment_tags']):
[371]147        v1 = seg[0]
148        v2 = seg[1]
149        for key in [(v1,v2),(v2,v1)]:
[969]150            if sides.has_key(key) and tag <> "":
[686]151                #"" represents null.  Don't put these into the dictionary
[969]152                #this creates a dict of lists of faces, indexed by tag
153                #tagged_edges.setdefault(tag,[]).append(sides[key])
154                tag_dict[sides[key]] = tag
[686]155
[969]156    return tag_dict
[371]157
[1556]158
[371]159def calc_sides(triangles):
160    #Build dictionary mapping from sides (2-tuple of points)
[1556]161    #to left hand side neighbouring triangle
[371]162    sides = {}
163    for id, triangle in enumerate(triangles):
164        a = triangle[0]
165        b = triangle[1]
[1556]166        c = triangle[2]
[371]167        sides[a,b] = (id, 2) #(id, face)
168        sides[b,c] = (id, 0) #(id, face)
169        sides[c,a] = (id, 1) #(id, face)
[1556]170    return sides
Note: See TracBrowser for help on using the repository browser.