Source code for meshpy.triangle

from meshpy.common import MeshInfoBase, dump_array
import meshpy._internals as internals

[docs]class MeshInfo(internals.TriMeshInfo, MeshInfoBase): _constituents = [ "points", "point_attributes", "point_markers", "elements", "element_attributes", "element_volumes", "neighbors", "facets", "facet_markers", "holes", "regions", "faces", "face_markers", "normals", ] def __getstate__(self): return self.number_of_point_attributes, \ self.number_of_element_attributes, \ [(name, getattr(self, name)) for name in self._constituents] def __setstate__(self, xxx_todo_changeme): (p_attr_count, e_attr_count, state) = xxx_todo_changeme self.number_of_point_attributes = p_attr_count self.number_of_element_attributes = e_attr_count for name, array in state: if name not in self._constituents: raise RuntimeError("Unknown constituent during unpickling") dest_array = getattr(self, name) if array is None: dest_array.deallocate() else: if len(dest_array) != len(array): dest_array.resize(len(array)) if not dest_array.allocated and len(array) > 0: dest_array.setup() for i, tup in enumerate(array): for j, v in enumerate(tup): dest_array[i, j] = v
[docs] def set_facets(self, facets, facet_markers=None): self.facets.resize(len(facets)) for i, facet in enumerate(facets): self.facets[i] = facet if facet_markers is not None: self.facet_markers.setup() for i, mark in enumerate(facet_markers): self.facet_markers[i] = mark
def dump(self): for name in self._constituents: dump_array(name, getattr(self, name))
[docs]def subdivide_facets(subdivisions, points, facets, facet_markers=None): """Return a new facets array in which the original facets are each subdivided into C{subdivisions} subfacets. This routine is useful if you have to prohibit the insertion of Steiner points on the boundary of your triangulation to allow the mesh to conform either to itself periodically or another given mesh. In this case, you may use this routine to create the necessary resolution along the boundary in a predefined way. @arg subdivisions: Either an C{int}, indicating a uniform number of subdivisions throughout, or a list of the same length as C{facets}, specifying a subdivision count for each individual facet. @arg points: A list of points referred to from the facets list. @arg facets: The list of old facets, in the form C{[(p1, p2), (p3,p4), ...]}. @arg facet_markers: Either C{None} or a list of facet markers of the same length as C{facets}. @return: The new tuple C{(new_points, new_facets)}. (Or C{(new_points, new_facets, new_facet_markers)} if C{facet_markers} is not C{None}.) """ def intermediate_points(pa, pb, n): for i in range(1, n): tau = i/n yield [pai*(1-tau) + tau*pbi for pai, pbi in zip(pa, pb)] if isinstance(subdivisions, int): from itertools import repeat subdiv_it = repeat(subdivisions, len(facets)) else: assert len(facets) == len(subdivisions) subdiv_it = subdivisions.__iter__() new_points = points[:] new_facets = [] if facet_markers is not None: assert len(facets) == len(facet_markers) new_facet_markers = [] for facet_idx, ((pidx_a, pidx_b), subdiv) in enumerate(zip(facets, subdiv_it)): facet_points = [pidx_a] for p in intermediate_points(points[pidx_a], points[pidx_b], subdiv): facet_points.append(len(new_points)) new_points.append(p) facet_points.append(pidx_b) for i, p1 in enumerate(facet_points[:-1]): p2 = facet_points[i+1] new_facets.append((p1, p2)) if facet_markers is not None: new_facet_markers.append(facet_markers[facet_idx]) if facet_markers is not None: return new_points, new_facets, new_facet_markers else: return new_points, new_facets
[docs]def build(mesh_info, verbose=False, refinement_func=None, attributes=False, volume_constraints=False, max_volume=None, allow_boundary_steiner=True, allow_volume_steiner=True, quality_meshing=True, generate_edges=None, generate_faces=False, min_angle=None, mesh_order=None, generate_neighbor_lists=False): """Triangulate the domain given in `mesh_info'.""" opts = "pzj" if quality_meshing: if min_angle is not None: opts += "q%f" % min_angle else: opts += "q" if mesh_order is not None: opts += "o%d" % mesh_order if verbose: opts += "VV" else: opts += "Q" if attributes: opts += "A" if volume_constraints: opts += "a" if max_volume: opts += "a%.20f" % max_volume if refinement_func is not None: opts += "u" if generate_edges is not None: from warnings import warn warn("generate_edges is deprecated--use generate_faces instead") generate_faces = generate_edges if generate_neighbor_lists is not None: opts += "n" if generate_faces: opts += "e" if not allow_volume_steiner: opts += "YY" if allow_boundary_steiner: raise ValueError("cannot allow boundary Steiner points when volume " "Steiner points are forbidden") else: if not allow_boundary_steiner: opts += "Y" # restore "C" locale--otherwise triangle might mis-parse stuff like "a0.01" try: import locale except ImportError: have_locale = False else: have_locale = True prev_num_locale = locale.getlocale(locale.LC_NUMERIC) locale.setlocale(locale.LC_NUMERIC, "C") try: mesh = MeshInfo() internals.triangulate(opts, mesh_info, mesh, MeshInfo(), refinement_func) finally: # restore previous locale if we've changed it if have_locale: locale.setlocale(locale.LC_NUMERIC, prev_num_locale) return mesh
[docs]def refine(input_p, verbose=False, refinement_func=None, quality_meshing=True, min_angle=None, generate_neighbor_lists=False): opts = "razj" if quality_meshing: if min_angle is not None: opts += "q%f" % min_angle else: opts += "q" if len(input_p.faces) != 0: opts += "p" if verbose: opts += "VV" else: opts += "Q" if refinement_func is not None: opts += "u" if generate_neighbor_lists is not None: opts += "n" output_p = MeshInfo() internals.triangulate(opts, input_p, output_p, MeshInfo(), refinement_func) return output_p
[docs]def write_gnuplot_mesh(filename, out_p, facets=False): gp_file = open(filename, "w") if facets: segments = out_p.facets else: segments = out_p.elements for points in segments: for pt in points: gp_file.write("%f %f\n" % tuple(out_p.points[pt])) gp_file.write("%f %f\n\n" % tuple(out_p.points[points[0]]))