From 97fae3184c70000a88ba66fb3798900d2c25466a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 26 Jul 2007 00:14:03 -0500 Subject: [PATCH] Usable periodicity in Tetgen. Fix index error reporting. Faces/facets cleanup. --- Makefile.in | 12 +- src/cpp/foreign_array_wrap.hpp | 18 +- src/cpp/wrap_tetgen.cpp | 363 +++++++++++++++++---------------- src/cpp/wrap_triangle.cpp | 36 ++-- src/python/common.py | 2 +- src/python/tet.py | 149 +++++++++++--- src/python/triangle.py | 20 +- test/test_tet_torus.py | 6 +- test/test_tetgen.py | 2 +- test/test_tetgen_2.py | 2 +- test/test_triangle.py | 2 +- 11 files changed, 365 insertions(+), 247 deletions(-) diff --git a/Makefile.in b/Makefile.in index 05f7fcc..f5761ac 100644 --- a/Makefile.in +++ b/Makefile.in @@ -1,8 +1,14 @@ -all: +.PHONY : all install clean tags + +all: tags ${PYTHON_EXE} setup.py build -install: +install: tags ${PYTHON_EXE} setup.py install clean: - rm build -Rf + rm -Rf build + rm tags + +tags: + ctags -R src || true diff --git a/src/cpp/foreign_array_wrap.hpp b/src/cpp/foreign_array_wrap.hpp index 56cd426..2a9ef4c 100644 --- a/src/cpp/foreign_array_wrap.hpp +++ b/src/cpp/foreign_array_wrap.hpp @@ -39,7 +39,8 @@ namespace { static object getitem(FA &self, long idx) { if (idx < 0) idx += self.size(); - if (idx >= (long) self.size()) PYTHON_ERROR(IndexError, "index out of bounds"); + if (idx < 0 || idx >= (long) self.size()) + PYTHON_ERROR(IndexError, "index out of bounds"); if (self.unit() > 1) { @@ -60,9 +61,11 @@ namespace { long i_sub = extract(idx[1]); if (i_main < 0) idx += self.size(); - if (i_main >= (long) self.size()) PYTHON_ERROR(IndexError, "index out of bounds"); + if (i_main < 0 || i_main >= (long) self.size()) + PYTHON_ERROR(IndexError, "index out of bounds"); if (i_sub < 0) idx += self.unit(); - if (i_sub >= (long) self.unit()) PYTHON_ERROR(IndexError, "subindex out of bounds"); + if (i_sub < 0 || i_sub >= (long) self.unit()) + PYTHON_ERROR(IndexError, "subindex out of bounds"); return object(self.getSub(i_main, i_sub)); } @@ -70,7 +73,8 @@ namespace { static void setitem(FA &self, long idx, object value) { if (idx < 0) idx += self.size(); - if (idx >= (long) self.size()) PYTHON_ERROR(IndexError, "index out of bounds"); + if (idx < 0 || idx >= (long) self.size()) + PYTHON_ERROR(IndexError, "index out of bounds"); if (self.unit() > 1) { @@ -92,9 +96,11 @@ namespace { long i_sub = extract(idx[1]); if (i_main < 0) idx += self.size(); - if (i_main >= (long) self.size()) PYTHON_ERROR(IndexError, "index out of bounds"); + if (i_main < 0 || i_main >= (long) self.size()) + PYTHON_ERROR(IndexError, "index out of bounds"); if (i_sub < 0) idx += self.unit(); - if (i_sub >= (long) self.unit()) PYTHON_ERROR(IndexError, "subindex out of bounds"); + if (i_main < 0 || i_sub >= (long) self.unit()) + PYTHON_ERROR(IndexError, "subindex out of bounds"); self.setSub(i_main, i_sub, v); } diff --git a/src/cpp/wrap_tetgen.cpp b/src/cpp/wrap_tetgen.cpp index 98c1fa4..b46b12a 100644 --- a/src/cpp/wrap_tetgen.cpp +++ b/src/cpp/wrap_tetgen.cpp @@ -15,238 +15,241 @@ using namespace std; -struct tMeshInfo : public tetgenio, public boost::noncopyable +namespace { - public: - tForeignArray Points; // in/out - tForeignArray PointAttributes; // in/out - tForeignArray PointMetricTensors; // in/out - tForeignArray PointMarkers; // in/out - - tForeignArray Elements; // in/out - tForeignArray ElementAttributes; // in/out - tForeignArray ElementVolumes; // out - tForeignArray Neighbors; // out - - tForeignArray Facets; - tForeignArray FacetMarkers; - - tForeignArray Holes; - tForeignArray Regions; - - tForeignArray FacetConstraints; - tForeignArray SegmentConstraints; + struct tMeshInfo : public tetgenio, public boost::noncopyable + { + public: + tForeignArray Points; // in/out + tForeignArray PointAttributes; // in/out + tForeignArray PointMetricTensors; // in/out + tForeignArray PointMarkers; // in/out + + tForeignArray Elements; // in/out + tForeignArray ElementAttributes; // in/out + tForeignArray ElementVolumes; // out + tForeignArray Neighbors; // out + + tForeignArray Facets; + tForeignArray FacetMarkers; + + tForeignArray Holes; + tForeignArray Regions; + + tForeignArray FacetConstraints; + tForeignArray SegmentConstraints; + + tForeignArray PBCGroups; + + tForeignArray Faces; + tForeignArray AdjacentElements; + tForeignArray FaceMarkers; + + tForeignArray Edges; + tForeignArray EdgeMarkers; + + public: + tMeshInfo() + : Points(pointlist, numberofpoints, 3), + PointAttributes(pointattributelist, numberofpoints, 0, &Points), + PointMetricTensors(pointmtrlist, numberofpoints, 0, &Points), + PointMarkers(pointmarkerlist, numberofpoints, 1, &Points), + + Elements(tetrahedronlist, numberoftetrahedra, 4), + ElementAttributes(tetrahedronattributelist, + numberoftetrahedra, 0, &Elements), + ElementVolumes(tetrahedronvolumelist, numberoftetrahedra, 1, &Elements), + Neighbors(neighborlist, numberoftetrahedra, 4, &Elements), + + Facets(facetlist, numberoffacets), + FacetMarkers(facetmarkerlist, numberoffacets, 1, &Facets), + + Holes(holelist, numberofholes, 3), + + Regions(regionlist, numberofregions, 5), + + FacetConstraints(facetconstraintlist, numberoffacetconstraints, 2), + SegmentConstraints(facetconstraintlist, numberofsegmentconstraints, 3), + + PBCGroups(pbcgrouplist, numberofpbcgroups), + + Faces(trifacelist, numberoftrifaces, 3), + AdjacentElements(adjtetlist, numberoftrifaces, 2, &Faces), + FaceMarkers(trifacemarkerlist, numberoftrifaces, 1, &Faces), + + Edges(edgelist, numberofedges, 2), + EdgeMarkers(edgemarkerlist, numberofedges, 1, &Edges) + { + } + + unsigned numberOfPointAttributes() const + { + return numberofpointattributes; + } + + unsigned numberOfPointMetricTensors() const + { + return numberofpointmtrs; + } + + unsigned numberOfElementAttributes() const + { + return numberoftetrahedronattributes; + } + + void setNumberOfPointAttributes(unsigned attrs) + { + PointAttributes.setUnit(attrs); + numberofpointattributes = attrs; + } + + void setNumberOfPointMetricTensors(unsigned mtrs) + { + PointMetricTensors.setUnit(mtrs); + numberofpointmtrs = mtrs; + } + + void setNumberOfElementAttributes(unsigned attrs) + { + ElementAttributes.setUnit(attrs); + numberoftetrahedronattributes = attrs; + } - tForeignArray PBCGroups; + /* + tTriangulationParameters &operator=(const tTriangulationParameters &src) + { + numberofpointattributes = src.numberofpointattributes ; + numberofcorners = src.numberofcorners; + numberoftriangleattributes = src.numberoftriangleattributes; - tForeignArray Faces; - tForeignArray AdjacentElements; - tForeignArray FaceMarkers; + Points = src.Points; + PointAttributes = src.PointAttributes; + PointMarkers = src.PointMarkers; - tForeignArray Edges; - tForeignArray EdgeMarkers; + Triangles = src.Triangles; + TriangleAttributes = src.TriangleAttributes; + TriangleAreas = src.TriangleAreas; + Neighbors = src.Neighbors; - public: - tMeshInfo() - : Points(pointlist, numberofpoints, 3), - PointAttributes(pointattributelist, numberofpoints, 0, &Points), - PointMetricTensors(pointmtrlist, numberofpoints, 0, &Points), - PointMarkers(pointmarkerlist, numberofpoints, 1, &Points), + Segments = src.Segments; + SegmentMarkers = src.SegmentMarkers; - Elements(tetrahedronlist, numberoftetrahedra, 4), - ElementAttributes(tetrahedronattributelist, - numberoftetrahedra, 0, &Elements), - ElementVolumes(tetrahedronvolumelist, numberoftetrahedra, 1, &Elements), - Neighbors(neighborlist, numberoftetrahedra, 4, &Elements), + Holes = src.Holes; - Facets(facetlist, numberoffacets), - FacetMarkers(facetmarkerlist, numberoffacets, 1, &Facets), + Regions = src.Regions; - Holes(holelist, numberofholes, 3), + Edges = src.Edges; + EdgeMarkers = src.EdgeMarkers; + Normals = src.Normals; - Regions(regionlist, numberofregions, 5), + return *this; + } + */ + }; - FacetConstraints(facetconstraintlist, numberoffacetconstraints, 2), - SegmentConstraints(facetconstraintlist, numberofsegmentconstraints, 3), - PBCGroups(pbcgrouplist, numberofpbcgroups), - Faces(trifacelist, numberoftrifaces, 3), - AdjacentElements(adjtetlist, numberoftrifaces, 2, &Faces), - FaceMarkers(trifacemarkerlist, numberoftrifaces, 1, &Faces), - Edges(edgelist, numberofedges, 2), - EdgeMarkers(edgemarkerlist, numberofedges, 1, &Edges) - { - } - - unsigned numberOfPointAttributes() const - { - return numberofpointattributes; - } + /* + tTriangulationParameters *copyTriangulationParameters(const tTriangulationParameters &src) + { + auto_ptr copy(new tTriangulationParameters); + *copy = src; + return copy.release(); + } + */ - unsigned numberOfPointMetricTensors() const - { - return numberofpointmtrs; - } - unsigned numberOfElementAttributes() const - { - return numberoftetrahedronattributes; - } - void setNumberOfPointAttributes(unsigned attrs) - { - PointAttributes.setUnit(attrs); - numberofpointattributes = attrs; - } - void setNumberOfPointMetricTensors(unsigned mtrs) + void tetrahedralizeWrapper(tetgenbehavior &bhv, tMeshInfo &in, tMeshInfo &out) + { + try { - PointMetricTensors.setUnit(mtrs); - numberofpointmtrs = mtrs; + tetrahedralize(&bhv, &in, &out); } - - void setNumberOfElementAttributes(unsigned attrs) + catch (int &i) { - ElementAttributes.setUnit(attrs); - numberoftetrahedronattributes = attrs; + throw runtime_error("TetGen runtime error code "+boost::lexical_cast(i)); } - /* - tTriangulationParameters &operator=(const tTriangulationParameters &src) - { - numberofpointattributes = src.numberofpointattributes ; - numberofcorners = src.numberofcorners; - numberoftriangleattributes = src.numberoftriangleattributes; - - Points = src.Points; - PointAttributes = src.PointAttributes; - PointMarkers = src.PointMarkers; - - Triangles = src.Triangles; - TriangleAttributes = src.TriangleAttributes; - TriangleAreas = src.TriangleAreas; - Neighbors = src.Neighbors; - - Segments = src.Segments; - SegmentMarkers = src.SegmentMarkers; - - Holes = src.Holes; - - Regions = src.Regions; - - Edges = src.Edges; - EdgeMarkers = src.EdgeMarkers; - Normals = src.Normals; - - return *this; - } - */ -}; + out.PointAttributes.fixUnit(out.numberofpointattributes); + out.PointMetricTensors.fixUnit(out.numberofpointmtrs); + out.ElementAttributes.fixUnit(out.numberoftetrahedronattributes); + } -/* -tTriangulationParameters *copyTriangulationParameters(const tTriangulationParameters &src) -{ - auto_ptr copy(new tTriangulationParameters); - *copy = src; - return copy.release(); -} -*/ + template + struct manage_new_internal_reference + : with_custodian_and_ward_postcall<0, owner_arg, Base> + { + typedef manage_new_object result_converter; + }; -void tetrahedralizeWrapper(tetgenbehavior &bhv, tMeshInfo &in, tMeshInfo &out) -{ - try + tForeignArray *facet_get_polygons(tetgenio::facet &self) { - tetrahedralize(&bhv, &in, &out); + return new tForeignArray( + self.polygonlist, self.numberofpolygons); } - catch (int &i) - { - throw runtime_error("TetGen runtime error code "+boost::lexical_cast(i)); - } - - out.PointAttributes.fixUnit(out.numberofpointattributes); - out.PointMetricTensors.fixUnit(out.numberofpointmtrs); - out.ElementAttributes.fixUnit(out.numberoftetrahedronattributes); -} -template -struct manage_new_internal_reference - : with_custodian_and_ward_postcall<0, owner_arg, Base> -{ - typedef manage_new_object result_converter; -}; - - - - -tForeignArray *facet_get_polygons(tetgenio::facet &self) -{ - return new tForeignArray( - self.polygonlist, self.numberofpolygons); -} - - - - -tForeignArray *facet_get_holes(tetgenio::facet &self) -{ - return new tForeignArray(self.holelist, self.numberofholes); -} + tForeignArray *facet_get_holes(tetgenio::facet &self) + { + return new tForeignArray(self.holelist, self.numberofholes); + } -tForeignArray *polygon_get_vertices(tetgenio::polygon &self) -{ - return new tForeignArray(self.vertexlist, self.numberofvertices); -} + tForeignArray *polygon_get_vertices(tetgenio::polygon &self) + { + return new tForeignArray(self.vertexlist, self.numberofvertices); + } -REAL pbcgroup_get_transmat_entry(tetgenio::pbcgroup &self, long i, long j) -{ - if (i < 0) i += 4; - if (j < 0) j += 4; + REAL pbcgroup_get_transmat_entry(tetgenio::pbcgroup &self, long i, long j) + { + if (i < 0) i += 4; + if (j < 0) j += 4; - if (i < 0 || i >= 4 || j < 0 || j >= 4) - PYTHON_ERROR(IndexError, "transform matrix index out of bounds"); - return self.transmat[i][j]; -} + if (i < 0 || i >= 4 || j < 0 || j >= 4) + PYTHON_ERROR(IndexError, "transform matrix index out of bounds"); + return self.transmat[i][j]; + } -void pbcgroup_set_transmat_entry(tetgenio::pbcgroup &self, long i, long j, REAL value) -{ - if (i < 0) i += 4; - if (j < 0) j += 4; + void pbcgroup_set_transmat_entry(tetgenio::pbcgroup &self, long i, long j, REAL value) + { + if (i < 0) i += 4; + if (j < 0) j += 4; - if (i < 0 || i >= 4 || j < 0 || j >= 4) - PYTHON_ERROR(IndexError, "transform matrix index out of bounds"); - self.transmat[i][j] = value; -} + if (i < 0 || i >= 4 || j < 0 || j >= 4) + PYTHON_ERROR(IndexError, "transform matrix index out of bounds"); + self.transmat[i][j] = value; + } -tForeignArray *pbcgroup_get_pointpairs(tetgenio::pbcgroup &self) -{ - return new tForeignArray(self.pointpairlist, self.numberofpointpairs); + tForeignArray *pbcgroup_get_pointpairs(tetgenio::pbcgroup &self) + { + return new tForeignArray(self.pointpairlist, self.numberofpointpairs, 2); + } } @@ -333,15 +336,18 @@ BOOST_PYTHON_MODULE(_tetgen) { typedef tetgenio::facet cl; class_("Facet", no_init) - .def("get_polygons", facet_get_polygons, manage_new_internal_reference<>()) - .def("get_holes", facet_get_holes, manage_new_internal_reference<>()) + .add_property("polygons", + make_function(facet_get_polygons, manage_new_internal_reference<>())) + .add_property("holes", + make_function(facet_get_holes, manage_new_internal_reference<>())) ; } { typedef tetgenio::polygon cl; class_("Polygon", no_init) - .def("get_vertices", polygon_get_vertices, manage_new_internal_reference<>()) + .add_property("vertices", + make_function(polygon_get_vertices, manage_new_internal_reference<>())) ; } @@ -352,7 +358,8 @@ BOOST_PYTHON_MODULE(_tetgen) .def_readwrite("facet_marker_2", &cl::fmark2) .def("get_transmat_entry", pbcgroup_get_transmat_entry) .def("set_transmat_entry", pbcgroup_set_transmat_entry) - .def("get_point_pairs", pbcgroup_get_pointpairs, manage_new_internal_reference<>()) + .add_property("point_pairs", + make_function(pbcgroup_get_pointpairs, manage_new_internal_reference<>())) ; } diff --git a/src/cpp/wrap_triangle.cpp b/src/cpp/wrap_triangle.cpp index 071f619..fa769d7 100644 --- a/src/cpp/wrap_triangle.cpp +++ b/src/cpp/wrap_triangle.cpp @@ -25,15 +25,15 @@ struct tMeshInfo : public triangulateio, public boost::noncopyable tForeignArray ElementVolumes; // in only tForeignArray Neighbors; // out only - tForeignArray Faces; // in/out - tForeignArray FaceMarkers; // in/out + tForeignArray Facets; // in/out + tForeignArray FacetMarkers; // in/out tForeignArray Holes; // in only tForeignArray Regions; // in only - tForeignArray Edges; // out only - tForeignArray EdgeMarkers; // out only + tForeignArray Faces; // out only + tForeignArray FaceMarkers; // out only tForeignArray Normals; // out only public: @@ -50,16 +50,16 @@ struct tMeshInfo : public triangulateio, public boost::noncopyable Neighbors(neighborlist, numberoftriangles, 3, &Elements, true), - Faces(segmentlist, numberofsegments, 2, NULL, true), - FaceMarkers(segmentmarkerlist, numberofsegments, 1, &Faces, true), + Facets(segmentlist, numberofsegments, 2, NULL, true), + FacetMarkers(segmentmarkerlist, numberofsegments, 1, &Facets, true), Holes(holelist, numberofholes, 2, NULL, true), Regions(regionlist, numberofregions, 4, NULL, true), - Edges(edgelist, numberofedges, 2, NULL, true), - EdgeMarkers(edgemarkerlist, numberofedges, 1, &Edges, true), - Normals(normlist, numberofedges, 2, &Edges, true) + Faces(edgelist, numberofedges, 2, NULL, true), + FaceMarkers(edgemarkerlist, numberofedges, 1, &Faces, true), + Normals(normlist, numberofedges, 2, &Faces, true) { numberofpointattributes = 0; numberofcorners = 3; @@ -103,15 +103,15 @@ struct tMeshInfo : public triangulateio, public boost::noncopyable ElementVolumes = src.ElementVolumes; Neighbors = src.Neighbors; - Faces = src.Faces; - FaceMarkers = src.FaceMarkers; + Facets = src.Facets; + FacetMarkers = src.FacetMarkers; Holes = src.Holes; Regions = src.Regions; - Edges = src.Edges; - EdgeMarkers = src.EdgeMarkers; + Faces = src.Faces; + FaceMarkers = src.FaceMarkers; Normals = src.Normals; return *this; @@ -205,15 +205,15 @@ BOOST_PYTHON_MODULE(_triangle) .def_readonly("element_volumes", &cl::ElementVolumes) .def_readonly("neighbors", &cl::Neighbors) - .def_readonly("faces", &cl::Faces) - .def_readonly("face_markers", &cl::FaceMarkers) + .def_readonly("facets", &cl::Facets) + .def_readonly("facet_markers", &cl::FacetMarkers) .def_readonly("holes", &cl::Holes) .def_readonly("regions", &cl::Regions) - .def_readonly("edges", &cl::Edges) - .def_readonly("edge_markers", &cl::EdgeMarkers) + .def_readonly("faces", &cl::Faces) + .def_readonly("face_markers", &cl::FaceMarkers) .def_readonly("normals", &cl::Normals) @@ -232,7 +232,7 @@ BOOST_PYTHON_MODULE(_triangle) exposePODForeignArray("RealArray"); exposePODForeignArray("IntArray"); - class_, tVertex, boost::noncopyable>("Vertex", no_init) + class_("Vertex", no_init) .add_property("x", &tVertex::x) .add_property("y", &tVertex::y) ; diff --git a/src/python/common.py b/src/python/common.py index 5bbe6f0..af74a01 100644 --- a/src/python/common.py +++ b/src/python/common.py @@ -160,7 +160,7 @@ class MeshInfoBase: frozenset([el[1], el[2], el[3]]), frozenset([el[2], el[0], el[3]]), ] - for fi, face in enumerate(tet_faces): + for fi, face in enumerate(faces): face2el.setdefault(face, []).append((ti, fi+1)) else: diff --git a/src/python/tet.py b/src/python/tet.py index 7a72d6d..c6fb238 100644 --- a/src/python/tet.py +++ b/src/python/tet.py @@ -6,7 +6,7 @@ import meshpy._tetgen as internals class MeshInfo(internals.MeshInfo, MeshInfoBase): - def set_faces(self, facets, markers=None): + def set_facets(self, facets, markers=None): if markers: assert len(markers) == len(facets) @@ -26,10 +26,23 @@ class MeshInfo(internals.MeshInfo, MeshInfoBase): for i, mark in enumerate(markers): self.facet_markers[i] = mark + @property + def face_vertex_indices_to_face_marker(self): + try: + return self._fvi2fm + except AttributeError: + result = {} + + for i, face in enumerate(self.faces): + result[frozenset(face)] = self.face_markers[i] + + self._fvi2fm = result + return result + def dump(self): for name in ["points"]: dump_array(name, getattr(self, name)) - for ifacet, facet in enumerate(self.facets): + for ifacet, facet in enumerate(self.faces): print "facet %d:" % ifacet for ipolygon, polygon in enumerate(facet.polygons): print " polygon %d: vertices [%s]" % \ @@ -59,29 +72,52 @@ class Options(internals.Options): def _PBCGroup_get_transmat(self): - return PBCGroupTransmat(self) + import pylinear.array as num + + return num.array( + [[self.get_transmat_entry(i,j) + for j in xrange(4)] + for i in xrange(4)]) + + +def _PBCGroup_set_transmat(self, matrix): + import pylinear.array as num + for i in xrange(4): + for j in xrange(4): + self.set_transmat_entry(i, j, matrix[i,j]) -class PBCGroupTransmat: - def __init__(self, pbcgroup): - self.pbcgroup = pbcgroup - def __getitem__(self, (i,j)): - return self.pcgroup.get_transmat_entry(i, j) - def __setitem__(self, (i,j), v): - return self.pcgroup.set_transmat_entry(i, j, v) +def _PBCGroup_set_transform(self, matrix=None, translation=None): + for i in xrange(4): + for j in xrange(4): + self.set_transmat_entry(i, j, 0) + self.set_transmat_entry(3, 3, 1) + if matrix is not None: + for i in xrange(3): + for j in xrange(3): + self.set_transmat_entry(i, j, matrix[i][j]) + else: + for i in xrange(3): + self.set_transmat_entry(i, i, 1) -internals.Facet.polygons = property(internals.Facet.get_polygons) -internals.Facet.holes = property(internals.Facet.get_holes) -internals.Polygon.vertices = property(internals.Polygon.get_vertices) -internals.PBCGroup.point_pairs = property(internals.PBCGroup.get_point_pairs) -internals.PBCGroup.matrix = property(_PBCGroup_get_transmat) + if translation is not None: + for i in xrange(3): + self.set_transmat_entry(i, 3, translation[i]) + + + + +internals.PBCGroup.matrix = property( + _PBCGroup_get_transmat, + _PBCGroup_set_transmat) +internals.PBCGroup.set_transform = _PBCGroup_set_transform @@ -114,7 +150,46 @@ EXT_CLOSED_IN_RZ = 2 -def generate_extrusion(rz_points, base_shape, closure=EXT_OPEN, point_idx_offset=0): +def generate_extrusion(rz_points, base_shape, closure=EXT_OPEN, + point_idx_offset=0, ring_tags=None, closure_tag=0): + """Extrude a given connected `base_shape' (a list of (x,y) points) + along the z axis. For each step in the extrusion, the base shape + is multiplied by a radius and shifted in the z direction. Radius + and z offset are given by `rz_points', which is a list of + (r, z) tuples. + + Returns (points, facets), where points is a list of points + and facets is a list of tuples of indices into that point list, + each tuple specifying a polygon. If point_idx_offset is not + zero, these indices start at this number. There may be a third + return value, see `facet_markers' below. + + The extrusion proceeds by generating quadrilaterals connecting each + ring. If any given radius in `rz_points' is 0, triangle fans are + produced instead of quads to provide non-degenerate closure. + + If `closure' is EXT_OPEN, no efforts are made to put end caps on the + extrusion. + + Specifying `closure' as EXT_CLOSE_IN_Z is (almost) + equivalent to adding points with zero radius at the beginning + at end of the rz_points list. The z coordinates are equal to + the existing first and last points in that list. (This case is + handled more efficiently by just generating big flat facets. + Since these facets are always flat, triangle fans are + unnecessary.) + + If `closure' is EXT_CLOSED_IN_RZ, then a torus-like structure + is assumed and the last ring is just connected to the first. + + If `ring_tags' is not None, it is an list of tags added to each + ring. There should be len(rz_points)-1 entries in this list. + If rings are added because of closure options, they receive the + tag `closure_tag'. If `facet_markers' is given, this function + returns (points, facets, tags), where tags is is a list containing + a tag for each generated facet. + """ + assert len(rz_points) > 0 def gen_ring(r, z): @@ -153,37 +228,61 @@ def generate_extrusion(rz_points, base_shape, closure=EXT_OPEN, point_idx_offset my_pairs = pair_with_successor(ring) return [(a, b, c, d) for ((a,b), (d,c)) in zip(prev_pairs, my_pairs)] + + def add_polygons(new_polys, tag): + polygons.extend(new_polys) + + if ring_tags is not None: + tags.extend(len(new_polys)*[tag]) + points = [] polygons = [] + tags = [] first_r, z = rz_points[0] first_ring = prev_ring = gen_ring(first_r, z) prev_r = first_r if closure == EXT_CLOSE_IN_Z: - polygons.append(first_ring) + add_polygons([tuple(first_ring)], closure_tag) - for r, z in rz_points[1:]: + for i, (r, z) in enumerate(rz_points[1:]): ring = gen_ring(r, z) - polygons.extend(connect_ring(r, ring, prev_r, prev_ring)) + + if ring_tags is not None: + ring_tag = ring_tags[i] + else: + ring_tag = None + + add_polygons( + connect_ring(r, ring, prev_r, prev_ring), + ring_tag) prev_ring = ring prev_r = r if closure == EXT_CLOSE_IN_Z: - polygons.append(prev_ring[::-1]) + add_polygons([tuple(prev_ring[::-1])], closure_tag) if closure == EXT_CLOSED_IN_RZ: - polygons.extend(connect_ring(first_r, first_ring, prev_r, prev_ring)) + add_polygons(connect_ring( + first_r, first_ring, prev_r, prev_ring), + closure_tag) - return points, polygons + if ring_tags is not None: + return points, polygons, tags + else: + return points, polygons -def generate_surface_of_revolution(rz_points, closure=EXT_OPEN, radial_subdiv=16, point_idx_offset=0): +def generate_surface_of_revolution(rz_points, + closure=EXT_OPEN, radial_subdiv=16, + point_idx_offset=0, ring_tags=None, + closure_tag=0): from math import sin, cos, pi dphi = 2*pi/radial_subdiv base_shape = [(cos(dphi*i), sin(dphi*i)) for i in range(radial_subdiv)] return generate_extrusion(rz_points, base_shape, closure=closure, - point_idx_offset=point_idx_offset) - + point_idx_offset=point_idx_offset, + ring_tags=ring_tags, closure_tag=closure_tag) diff --git a/src/python/triangle.py b/src/python/triangle.py index 3a1792e..95c2c5c 100644 --- a/src/python/triangle.py +++ b/src/python/triangle.py @@ -9,10 +9,10 @@ class MeshInfo(internals.MeshInfo, MeshInfoBase): "points", "point_attributes", "point_markers", "elements", "element_attributes", "element_volumes", "neighbors", - "faces", "face_markers", + "facets", "facet_markers", "holes", "regions", - "edges", "edge_markers", + "faces", "face_markers", "normals", ] @@ -51,16 +51,16 @@ class MeshInfo(internals.MeshInfo, MeshInfoBase): for j,v in enumerate(tup): dest_array[i, j] = v - def set_faces(self, faces, face_markers=None): - self.faces.resize(len(faces)) + def set_facets(self, facets, facet_markers=None): + self.facets.resize(len(facets)) - for i, face in enumerate(faces): - self.faces[i] = face + for i, facet in enumerate(facets): + self.facets[i] = facet - if face_markers is not None: - self.face_markers.setup() - for i, mark in enumerate(face_markers): - self.face_markers[i] = mark + 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: diff --git a/test/test_tet_torus.py b/test/test_tet_torus.py index 9a1b504..6f1aea3 100644 --- a/test/test_tet_torus.py +++ b/test/test_tet_torus.py @@ -4,9 +4,9 @@ def main(): EXT_CLOSED_IN_RZ big_r = 3 - little_r = 1 + little_r = 2.9 - points = 30 + points = 50 dphi = 2*pi/points rz = [(big_r+little_r*cos(i*dphi), little_r*sin(i*dphi)) @@ -17,7 +17,7 @@ def main(): closure=EXT_CLOSED_IN_RZ, radial_subdiv=20) mesh_info.set_points(points) - mesh_info.set_facets(facets, [1 for i in range(len(facets))]) + mesh_info.set_facets(facets, len(facets) *[1]) mesh_info.save_nodes("torus") mesh_info.save_poly("torus") mesh = build(mesh_info) diff --git a/test/test_tetgen.py b/test/test_tetgen.py index 6cf2e89..27c8628 100644 --- a/test/test_tetgen.py +++ b/test/test_tetgen.py @@ -17,7 +17,7 @@ def main(): (0,2,12), ]) - mesh_info.set_faces([ + mesh_info.set_facets([ [0,1,2,3], [4,5,6,7], [0,4,5,1], diff --git a/test/test_tetgen_2.py b/test/test_tetgen_2.py index a51830b..7cd10f9 100644 --- a/test/test_tetgen_2.py +++ b/test/test_tetgen_2.py @@ -14,7 +14,7 @@ def main(): points, facets = generate_surface_of_revolution(simple_rz) mesh_info.set_points(points) - mesh_info.set_faces(facets, [0 for i in range(len(facets))]) + mesh_info.set_facets(facets, [0 for i in range(len(facets))]) #mesh_info.save_nodes("test") #mesh_info.save_poly("test") #mesh_info.load_poly("test") diff --git a/test/test_triangle.py b/test/test_triangle.py index 297eaed..6066de1 100644 --- a/test/test_triangle.py +++ b/test/test_triangle.py @@ -28,7 +28,7 @@ def needs_refinement( vert_origin, vert_destination, vert_apex, area ): info = triangle.MeshInfo() info.set_points(points) info.set_holes([(0,0)]) -info.set_faces(round_trip_connect(0, len(points)-1)) +info.set_facets(round_trip_connect(0, len(points)-1)) mesh = triangle.build(info, refinement_func=needs_refinement) -- GitLab