Package meshpy :: Module triangle
[hide private]
[frames] | no frames]

Source Code for Module meshpy.triangle

  1  from __future__ import division 
  2  from meshpy.common import MeshInfoBase, dump_array 
  3  import meshpy._triangle as internals 
  4   
  5   
  6   
  7   
8 -class MeshInfo(internals.MeshInfo, MeshInfoBase):
9 _constituents = [ 10 "points", "point_attributes", "point_markers", 11 "elements", "element_attributes", "element_volumes", 12 "neighbors", 13 "facets", "facet_markers", 14 "holes", 15 "regions", 16 "faces", "face_markers", 17 "normals", 18 ] 19
20 - def __getstate__(self):
21 def dump_array(array): 22 try: 23 return [ 24 [array[i,j] for j in range(array.unit)] 25 for i in range(len(array))] 26 except RuntimeError: 27 # not allocated 28 return None
29 30 return self.number_of_point_attributes, \ 31 self.number_of_element_attributes, \ 32 [(name, dump_array(getattr(self, name))) for name in self._constituents]
33
34 - def __setstate__(self, (p_attr_count, e_attr_count, state)):
35 self.number_of_point_attributes = p_attr_count 36 self.number_of_element_attributes = e_attr_count 37 for name, array in state: 38 if name not in self._constituents: 39 raise RuntimeError, "Unknown constituent during unpickling" 40 41 dest_array = getattr(self, name) 42 43 if array is None: 44 dest_array.deallocate() 45 else: 46 if len(dest_array) != len(array): 47 dest_array.resize(len(array)) 48 if not dest_array.allocated and len(array)>0: 49 dest_array.setup() 50 51 for i,tup in enumerate(array): 52 for j,v in enumerate(tup): 53 dest_array[i, j] = v
54
55 - def set_facets(self, facets, facet_markers=None):
56 self.facets.resize(len(facets)) 57 58 for i, facet in enumerate(facets): 59 self.facets[i] = facet 60 61 if facet_markers is not None: 62 self.facet_markers.setup() 63 for i, mark in enumerate(facet_markers): 64 self.facet_markers[i] = mark
65
66 - def dump(self):
67 for name in self._constituents: 68 dump_array(name, getattr(self, name))
69 70 71 72
73 -def subdivide_facets(subdivisions, points, facets, facet_markers=None):
74 """Return a new facets array in which the original facets are 75 each subdivided into C{subdivisions} subfacets. 76 77 This routine is useful if you have to prohibit the insertion of Steiner 78 points on the boundary of your triangulation to allow the mesh to conform 79 either to itself periodically or another given mesh. In this case, you may 80 use this routine to create the necessary resolution along the boundary 81 in a predefined way. 82 83 @arg subdivisions: Either an C{int}, indicating a uniform number of subdivisions 84 throughout, or a list of the same length as C{facets}, specifying a subdivision 85 count for each individual facet. 86 @arg points: A list of points referred to from the facets list. 87 @arg facets: The list of old facets, in the form C{[(p1, p2), (p3,p4), ...]}. 88 @arg facet_markers: Either C{None} or a list of facet markers of the same length 89 as C{facets}. 90 @return: The new tuple C{(new_points, new_facets)}. 91 (Or C{(new_points, new_facets, new_facet_markers)} if C{facet_markers} is not 92 C{None}.) 93 """ 94 95 def intermediate_points(pa, pb, n): 96 for i in range(1, n): 97 tau = i/n 98 yield [pai*(1-tau) + tau*pbi for pai, pbi in zip(pa, pb)]
99 100 if isinstance(subdivisions, int): 101 from itertools import repeat 102 subdiv_it = repeat(subdivisions, len(facets)) 103 else: 104 assert len(facets) == len(subdivisions) 105 subdiv_it = subdivisions.__iter__() 106 107 new_points = points[:] 108 new_facets = [] 109 110 if facet_markers is not None: 111 assert len(facets) == len(facet_markers) 112 new_facet_markers = [] 113 114 for facet_idx, ((pidx_a, pidx_b), subdiv) in enumerate(zip(facets, subdiv_it)): 115 facet_points = [pidx_a] 116 for p in intermediate_points(points[pidx_a], points[pidx_b], subdiv): 117 facet_points.append(len(new_points)) 118 new_points.append(p) 119 facet_points.append(pidx_b) 120 121 for i, p1 in enumerate(facet_points[:-1]): 122 p2 = facet_points[i+1] 123 new_facets.append((p1, p2)) 124 125 if facet_markers is not None: 126 new_facet_markers.append(facet_markers[facet_idx]) 127 128 if facet_markers is not None: 129 return new_points, new_facets, new_facet_markers 130 else: 131 return new_points, new_facets 132 133 134 135
136 -def build(mesh_info, verbose=False, refinement_func=None, attributes=False, 137 volume_constraints=False, max_volume=None, allow_boundary_steiner=True, 138 allow_volume_steiner=True, quality_meshing=True, 139 generate_edges=None, generate_faces=False, min_angle=None):
140 """Triangulate the domain given in `mesh_info'.""" 141 opts = "pzj" 142 if quality_meshing: 143 if min_angle is not None: 144 opts += "q%f" % min_angle 145 else: 146 opts += "q" 147 148 if verbose: 149 opts += "VV" 150 else: 151 opts += "Q" 152 153 if attributes: 154 opts += "A" 155 156 if volume_constraints: 157 opts += "a" 158 if max_volume: 159 raise ValueError, "cannot specify both volume_constraints and max_area" 160 elif max_volume: 161 opts += "a%s" % repr(max_volume) 162 163 if refinement_func is not None: 164 opts += "u" 165 166 if generate_edges is not None: 167 from warnings import warn 168 warn("generate_edges is deprecated--use generate_faces instead") 169 generate_faces = generate_edges 170 171 if generate_faces: 172 opts += "e" 173 174 if not allow_volume_steiner: 175 opts += "YY" 176 if allow_boundary_steiner: 177 raise ValueError("cannot allow boundary Steiner points when volume " 178 "Steiner points are forbidden") 179 else: 180 if not allow_boundary_steiner: 181 opts += "Y" 182 183 # restore "C" locale--otherwise triangle might mis-parse stuff like "a0.01" 184 try: 185 import locale 186 except ImportErorr: 187 have_locale = False 188 else: 189 have_locale = True 190 prev_num_locale = locale.getlocale(locale.LC_NUMERIC) 191 locale.setlocale(locale.LC_NUMERIC, "C") 192 193 try: 194 mesh = MeshInfo() 195 internals.triangulate(opts, mesh_info, mesh, MeshInfo(), refinement_func) 196 finally: 197 # restore previous locale if we've changed it 198 if have_locale: 199 locale.setlocale(locale.LC_NUMERIC, prev_num_locale) 200 201 return mesh
202 203 204 205
206 -def refine(input_p, verbose=False, refinement_func=None):
207 opts = "razj" 208 if len(input_p.faces) != 0: 209 opts += "p" 210 if verbose: 211 opts += "VV" 212 else: 213 opts += "Q" 214 if refinement_func is not None: 215 opts += "u" 216 output_p = MeshInfo() 217 internals.triangulate(opts, input_p, output_p, MeshInfo(), refinement_func) 218 return output_p
219 220 221 222
223 -def write_gnuplot_mesh(filename, out_p, facets=False):
224 gp_file = open(filename, "w") 225 226 if facets: 227 segments = out_p.facets 228 else: 229 segments = out_p.elements 230 231 for points in segments: 232 for pt in points: 233 gp_file.write("%f %f\n" % tuple(out_p.points[pt])) 234 gp_file.write("%f %f\n\n" % tuple(out_p.points[points[0]]))
235