[1] | 1 | def pmesh_to_domain(fileName, tags = None, setting_function = None): |
---|
| 2 | """ |
---|
| 3 | """ |
---|
| 4 | #FIXME: The current design doesn't seem to accomadate tags and |
---|
| 5 | # setting_functions being passed into the domain at this point. |
---|
| 6 | |
---|
| 7 | #FIXME: Plus the names of the functions are no longer appropriate, |
---|
| 8 | #since domain objects aren't being returned. |
---|
| 9 | |
---|
| 10 | import sys |
---|
| 11 | from config import pmesh_filename |
---|
| 12 | sys.path.append(pmesh_filename) |
---|
| 13 | from load_mesh.loadASCII import import_trianglulation |
---|
| 14 | |
---|
| 15 | try: |
---|
| 16 | meshdic = import_trianglulation(fileName) |
---|
| 17 | except IOError, e: |
---|
| 18 | msg = 'Could not open file %s ' %fileName |
---|
| 19 | raise IOError, msg |
---|
| 20 | |
---|
| 21 | |
---|
| 22 | return pmesh_dictionary_to_domain(meshdic, |
---|
| 23 | setting_function = setting_function) |
---|
| 24 | |
---|
| 25 | |
---|
| 26 | |
---|
| 27 | def pmesh_dictionary_to_domain(meshdic, setting_function = None): |
---|
| 28 | """ |
---|
| 29 | convert a pmesh dictionary to a list of Volumes. |
---|
| 30 | Also, return a list of triangles which have boundary tags |
---|
| 31 | meshdic structure; |
---|
| 32 | generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) |
---|
| 33 | generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...] |
---|
| 34 | generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers) |
---|
| 35 | generated segment marker list: [S1Marker, S2Marker, ...] (list of ints) |
---|
| 36 | triangle list: [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers) |
---|
| 37 | triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor |
---|
| 38 | triangle attribute list: [T1att, T2att, ...] (list of floats) |
---|
| 39 | """ |
---|
| 40 | |
---|
| 41 | vertex_coordinates = meshdic['generatedpointlist'] |
---|
| 42 | volumes = meshdic['generatedtrianglelist'] |
---|
| 43 | |
---|
| 44 | #if setting_function: |
---|
| 45 | # if not type(setting_function) is ListType: |
---|
| 46 | # setting_function = [setting_function] |
---|
| 47 | # for funct in setting_function: |
---|
| 48 | # meshdic = funct(meshdic, vertices = mesh_vertices, |
---|
| 49 | # volumes = volumes) |
---|
| 50 | |
---|
| 51 | marker_dict = pmesh_dict_to_marker_dict(volumes, meshdic, |
---|
| 52 | vertex_coordinates) |
---|
| 53 | |
---|
| 54 | return vertex_coordinates, volumes, marker_dict |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | #FIXME: The issue is whether this format should be stored in the tsh file |
---|
| 60 | #instead of having to be created here? |
---|
| 61 | |
---|
| 62 | #FIXME: Another issue is that the tsh file stores consecutive |
---|
| 63 | #indices explicitly. This is really redundant. |
---|
| 64 | #Suggest looking at obj and our own sww format and also consider |
---|
| 65 | #using netCDF. |
---|
| 66 | |
---|
| 67 | |
---|
| 68 | def pmesh_dict_to_marker_dict(triangles,meshdic,Vertices): |
---|
| 69 | """ Convert the pmesh dictionary (meshdic) description of boundary tags |
---|
| 70 | to a dictionary of markers, indexed with volume id and face number. |
---|
| 71 | """ |
---|
| 72 | sides = calc_sides(triangles) |
---|
| 73 | marker_dict = {} |
---|
| 74 | for seg, marker in map(None,meshdic['generatedsegmentlist'], |
---|
| 75 | meshdic['generatedsegmentmarkerlist']): |
---|
| 76 | v1 = seg[0] |
---|
| 77 | v2 = seg[1] |
---|
| 78 | for key in [(v1,v2),(v2,v1)]: |
---|
| 79 | if sides.has_key(key): |
---|
| 80 | #this creates a dict of lists of faces, indexed by marker |
---|
| 81 | #tagged_edges.setdefault(marker,[]).append(sides[key]) |
---|
| 82 | marker_dict[sides[key]] = marker |
---|
| 83 | return marker_dict |
---|
| 84 | |
---|
| 85 | |
---|
| 86 | def calc_sides(triangles): |
---|
| 87 | #Build dictionary mapping from sides (2-tuple of points) |
---|
| 88 | #to left hand side neighbouring triangle |
---|
| 89 | sides = {} |
---|
| 90 | for id, triangle in enumerate(triangles): |
---|
| 91 | a = triangle[0] |
---|
| 92 | b = triangle[1] |
---|
| 93 | c = triangle[2] |
---|
| 94 | sides[a,b] = (id, 2) #(id, face) |
---|
| 95 | sides[b,c] = (id, 0) #(id, face) |
---|
| 96 | sides[c,a] = (id, 1) #(id, face) |
---|
| 97 | return sides |
---|
| 98 | |
---|