From d62217060c811e0b5f0386164e595550ec5bdbe5 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 12 Feb 2007 13:29:35 -0500 Subject: [PATCH] [meshpy @ inform@tiker.net-20070212182935-55dbbd6e2f76c2c6] Fix FPU control word setting for Intel architectures. More identifier renaming fallout. Read Superfish format, generate Surfaces of Revolution. --- .bzrignore | 11 ++- meshpy/tet.py | 13 +++- meshpy/triangle.py | 4 +- setup.py | 1 + src/predicates.cpp | 8 +- src/triangle.c | 4 + src/wrap_tetgen.cpp | 13 ++++ test/gun.am | 73 ++++++++++++++++++ test/test_tetgen.py | 9 ++- test/test_tetgen_2.py | 174 ++++++++++++++++++++++++++++++++++++++++++ 10 files changed, 299 insertions(+), 11 deletions(-) create mode 100644 test/gun.am create mode 100644 test/test_tetgen_2.py diff --git a/.bzrignore b/.bzrignore index ed9a0f4..d5e7edc 100644 --- a/.bzrignore +++ b/.bzrignore @@ -1,4 +1,11 @@ build -*.node -*.poly *.dat +test/*.lua +MANIFEST +*.vtk +*.ele +*.face +*.poly +*.node +test/ParaViewTrace*.pvs +dist diff --git a/meshpy/tet.py b/meshpy/tet.py index 054b5d1..c8bb3db 100644 --- a/meshpy/tet.py +++ b/meshpy/tet.py @@ -34,6 +34,17 @@ class MeshInfo(internals.MeshInfo, MeshInfoBase): print " polygon %d: vertices [%s]" % \ (ipolygon, ",".join(str(vi) for vi in polygon.vertices)) + def write_vtk(self, filename): + import pyvtk + vtkelements = pyvtk.VtkData( + pyvtk.UnstructuredGrid( + self.points, + tetra=self.elements), + "Mesh") + vtkelements.tofile(filename) + + + class Options(internals.Options): @@ -62,5 +73,3 @@ def build(mesh_info, options=Options()): - - diff --git a/meshpy/triangle.py b/meshpy/triangle.py index 7247e78..80fc30b 100644 --- a/meshpy/triangle.py +++ b/meshpy/triangle.py @@ -58,7 +58,7 @@ class MeshInfo(internals.MeshInfo, MeshInfoBase): if segment_markers is not None: assert len(segment_markers) == len(segments) for i, mark in enumerate(segment_markers): - self.SegmentMarkers[i] = mark + self.segment_markers[i] = mark def dump(self): for name in self._constituents: @@ -88,7 +88,7 @@ def build(mesh_info, verbose=False, refinement_func=None): def refine(input_p, verbose=False, refinement_func=None): opts = "razj" - if input_p.Segments.size() != 0: + if len(input_p.segments) != 0: opts += "p" if verbose: opts += "VV" diff --git a/setup.py b/setup.py index e2ad3b3..2533c5c 100644 --- a/setup.py +++ b/setup.py @@ -17,6 +17,7 @@ triangle_macros = [ tetgen_macros = [ ("TETLIBRARY", 1), + ( "SELF_CHECK", 1 ) , ] setup(name="MeshPy", diff --git a/src/predicates.cpp b/src/predicates.cpp index 0d41e5b..ea9bf49 100644 --- a/src/predicates.cpp +++ b/src/predicates.cpp @@ -113,6 +113,10 @@ /* */ /*****************************************************************************/ +#if defined(__linux__) && defined(__i386__) + #define LINUX 1 +#endif + #include #include #include @@ -149,8 +153,8 @@ /* which is disastrously slow. A faster way on IEEE machines might be to */ /* mask the appropriate bit, but that's difficult to do in C. */ -#define Absolute(a) ((a) >= 0.0 ? (a) : -(a)) -/* #define Absolute(a) fabs(a) */ +/* #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))*/ +#define Absolute(a) fabs(a) /* Many of the operations are broken up into two pieces, a main part that */ /* performs an approximate operation, and a "tail" that computes the */ diff --git a/src/triangle.c b/src/triangle.c index eba2973..8ea99e7 100644 --- a/src/triangle.c +++ b/src/triangle.c @@ -269,6 +269,10 @@ /* #define CPU86 */ /* #define LINUX */ +#if defined(__linux__) && defined(__i386__) + #define LINUX 1 +#endif + #define INEXACT /* Nothing */ /* #define INEXACT volatile */ diff --git a/src/wrap_tetgen.cpp b/src/wrap_tetgen.cpp index 2d1dfbc..21b94d8 100644 --- a/src/wrap_tetgen.cpp +++ b/src/wrap_tetgen.cpp @@ -251,6 +251,19 @@ BOOST_PYTHON_MODULE(_tetgen) .DEF_METHOD(save_neighbors) .DEF_METHOD(save_poly) + .DEF_METHOD(load_node) + .DEF_METHOD(load_addnodes) + .DEF_METHOD(load_pbc) + .DEF_METHOD(load_var) + .DEF_METHOD(load_mtr) + .DEF_METHOD(load_poly) + .DEF_METHOD(load_off) + .DEF_METHOD(load_ply) + .DEF_METHOD(load_stl) + .DEF_METHOD(load_medit) + .DEF_METHOD(load_plc) + .DEF_METHOD(load_tetmesh) + /* .def("copy", ©TriangulationParameters, return_value_policy()) diff --git a/test/gun.am b/test/gun.am new file mode 100644 index 0000000..33cb4be --- /dev/null +++ b/test/gun.am @@ -0,0 +1,73 @@ +Ballistic bunch compression gun, 2.6-cell model, thermionic-type cells +!File received from John Lewellen +® kprob=1, freq= 2856 , +!xmax= 15 , ymax= 4.2 , +kmax = 300, lmax = 100, + +! drive point for cathode cell +xdri=.16, ydri=4.191 + +! drive point for full cell #1 +!xdri=5.848, ydri=4.191 + +! drive point for full cell #2 +!xdri=11.096, ydri=4.191 +& + + +! define the cathode cell +&po x=0., y=0. & +&po x= 0 , y= .6 & +&po x= .16 , y= 1.35 & +&po x=.16 , y=4.191 & +&po x=.6 , y=4.191 & +&po nt=2, x0= .6 , y0= 1.727 , r= 2.464 , theta=0. & +&po x=3.064 , y=1.4818 & +&po nt=2, x0=2.861 , y0=1.4818, r=.203, theta=270 & +&po x=2.75 , y=1.2788 & +&po nt=2, x0=2.750 , y0=1.0588, r=.22, theta=180 & +&po x=2.530 , y=1.008 & +&po nt=2, x0=3.03 , y0=1.008, r=.5, theta=270 & +&po x=3.224 , y=.508 & + + +! full cell 1 +&po x=3.481 , y=0.508 & +&po nt=2, x0=3.481 , y0=1.008 , r=0.5 , theta=0. & +&po x=3.981 , y=1.0588 & +&po nt=2, x0=3.761 , y0=1.0588 , r=0.22 , theta=90. & +&po x=3.587 , y=1.2788 & +&po nt=2, x0=3.587 , y0=1.4818 , r=0.203 , theta=180. & +&po x=3.384 , y=1.727 & +&po nt=2, x0=5.848 , y0=1.727 , r=2.464 , theta=90. & +&po nt=2, x0=5.858 , y0=1.727 , r=2.464 , theta= 0. & +&po x=8.312 , y=1.4818 & +&po nt=2, x0=8.109 , y0=1.4818 , r=0.203 , theta=270. & +&po x=7.935 , y=1.2788 & +&po nt=2, x0=7.935 , y0=1.0588 , r=0.22 , theta=180. & +&po x=7.715 , y=1.008 & +&po nt=2, x0=8.215 , y0=1.008 , r=0.5 , theta=270. & +&po x=8.472 , y=0.508 & + + + +! full cell 2 +&po x=8.729 , y=0.508 & +&po nt=2, x0=8.729 , y0=1.008 , r=0.5 , theta=0. & +&po x=9.229 , y=1.0588 & +&po nt=2, x0=9.009 , y0=1.0588 , r=0.22 , theta=90. & +&po x=8.835 , y=1.2788 & +&po nt=2, x0=8.835 , y0=1.4818 , r=0.203 , theta=180. & +&po x=8.632 , y=1.727 & +&po nt=2, x0=11.096 , y0=1.727 , r=2.464 , theta=90. & +&po nt=2, x0=11.096 , y0=1.727 , r=2.464 , theta= 0. & +&po x=13.560 , y=1.4818 & +&po nt=2, x0=13.357 , y0=1.4818 , r=0.203 , theta=270. & +&po x=13.183 , y=1.2788 & +&po nt=2, x0=13.183 , y0=1.0588 , r=0.22 , theta=180. & +&po x=12.963 , y=1.008 & +&po nt=2, x0=13.463 , y0=1.008 , r=0.5 , theta=270. & +&po x=13.72 , y=0.508 & +&po x=15.0 , y=0.508 & +&po x=15.0 , y=0. & +&po x=0, y=0 & diff --git a/test/test_tetgen.py b/test/test_tetgen.py index 005ba55..e7115a7 100644 --- a/test/test_tetgen.py +++ b/test/test_tetgen.py @@ -22,12 +22,15 @@ def main(): [1,5,6,2], [2,6,7,3], [3,7,4,0], - ], - [-1,-2,0,0,0,0]) + ]) #mesh_info.dump() #mesh_info.save_nodes("test") #mesh_info.save_poly("test") - build(mesh_info) + mesh = build(mesh_info) + #for el in mesh.elements: + #print el + mesh.write_vtk("test.vtk") + diff --git a/test/test_tetgen_2.py b/test/test_tetgen_2.py new file mode 100644 index 0000000..ac461b5 --- /dev/null +++ b/test/test_tetgen_2.py @@ -0,0 +1,174 @@ +from meshpy.tet import MeshInfo, build + + + + +def parse_superfish_format(filename): + import re + + lines = [line.strip() for line in file(filename, "r").readlines() + if line.strip().startswith("&po")] + + # parse the file + line_re = re.compile(r"^\&po\s+(.*)\s*\&$") + key_val_re = re.compile(r"^\s*([0-9a-zA-Z]+)\s*\=\s*(.*)\s*$") + proplines = [] + for line in lines: + line_match = line_re.match(line) + assert line_match + data = line_match.group(1).split(",") + + properties = {} + for datum in data: + datum_match = key_val_re.match(datum.strip()) + assert datum_match + properties[datum_match.group(1)] = float(datum_match.group(2)) + proplines.append(properties) + + # concoct x-y-points + from math import atan2, sin, cos, pi, ceil + import pylinear.array as num + import pylinear.computation as comp + + def add_point(pt): + points.append((pt[1], pt[0])) + + points = [] + for i, p in enumerate(proplines): + if "nt" in p: + # draw arc + last_point = points[-1] + center = num.array((p["x0"], p["y0"])) + r = p["r"] + theta = pi/180.*p["theta"] + + tangent = num.array((cos(theta), sin(theta))) + upward_normal = num.array((-sin(theta), cos(theta))) + + #print 180*theta/pi, last_point, center, tangent, upward_normal, last_point-center + if (last_point - center)*upward_normal < 0: + # draw ccw (positive) arc / upward + phi_start = 3*pi/2 + theta + phi_end = phi_start + pi/2 + else: + # draw cw (negative) arc / downward + phi_start = pi/2 + theta + phi_end = phi_start - pi/2 + steps = 10 + dphi = (phi_end-phi_start) / steps + phi = phi_start+dphi + + for i in range(steps): + points.append(center + r*num.array((cos(phi), sin(phi)))) + phi += dphi + else: + points.append(num.array((p["x"], p["y"]))) + + #outf = file("gun-lines.dat", "w") + #for p in points: + #outf.write("%f\t%f\n" % (p[0], p[1])) + #outf.close() + + # turn (x,y) into (r,z) + return [(p[1], p[0]) for p in points] + + + + + + + + +def generate_surface_of_revolution(rz_points, radial_subdiv=16, point_offset=0): + assert len(rz_points) > 0 + + from math import sin, cos, pi + + def gen_point(r, phi, z): + return (r*cos(phi), r*sin(phi), z) + + def gen_ring(r, z): + if r == 0: + p_indices = [p0+len(points)] + points.append(gen_point(r, 0, z)) + else: + p_indices = [p0+len(points)+i for i in range(radial_subdiv)] + points.extend([gen_point(r, dphi*i, z) for i in range(radial_subdiv)]) + return p_indices + + def pair_with_successor(l): + n = len(l) + return [(l[i], l[(i+1)%n]) for i in range(n)] + + p0 = point_offset + points = [] + polygons = [] + + dphi = 2*pi/radial_subdiv + + last_r, last_z = rz_points[0] + last_ring = gen_ring(last_r, last_z) + + for r, z in rz_points[1:]: + ring = gen_ring(r, z) + if last_r == 0: + # make opening fan + assert len(last_ring) == 1 + start_pt = last_ring[0] + if r != 0: + polygons.extend( + [(start_pt, succ, pt) for pt, succ in pair_with_successor(ring)] + ) + elif r == 0: + # make closing fan + assert len(ring) == 1 + end_pt = ring[0] + polygons.extend( + [(pt, succ, end_pt) for pt, succ in pair_with_successor(last_ring)] + ) + else: + # make quad strip + last_pairs = pair_with_successor(last_ring) + my_pairs = pair_with_successor(ring) + polygons.extend( + [(a, b, c, d) for ((a,b), (d,c)) in zip(last_pairs, my_pairs)] + ) + + last_ring = ring + last_r = r + last_z = z + + return points, polygons + + + + + + + +def main(): + simple_rz = [ + (0,0), + (1,1), + (1,2), + (0,3), + ] + mesh_info = MeshInfo() + rz_points = parse_superfish_format("gun.am") + points, facets = generate_surface_of_revolution(rz_points) + + mesh_info.set_points(points) + 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") + mesh = build(mesh_info) + mesh.write_vtk("my_mesh.vtk") + #mesh.save_elements("gun") + #mesh.save_nodes("gun") + + + + +if __name__ == "__main__": + main() -- GitLab