def test_tetgen():
    from meshpy.tet import MeshInfo, build
    mesh_info = MeshInfo()

    mesh_info.set_points([
        (0, 0, 0),
        (2, 0, 0),
        (2, 2, 0),
        (0, 2, 0),
        (0, 0, 12),
        (2, 0, 12),
        (2, 2, 12),
        (0, 2, 12),
        ])

    mesh_info.set_facets([
        [0, 1, 2, 3],
        [4, 5, 6, 7],
        [0, 4, 5, 1],
        [1, 5, 6, 2],
        [2, 6, 7, 3],
        [3, 7, 4, 0],
        ])

    build(mesh_info)

def main():
    #Test1: mesh testing
    MeshPoints = []
    MeshFacets = []
    with open('MeshFacets', 'rb') as input:
        MeshFacets = pickle.load(input)
    with open('MeshPoints', 'rb') as input:
        MeshPoints = pickle.load(input)

    mesh_info = MeshInfo()

    mesh_info.regions.resize(1)
    mesh_info.regions[0] = [
                            MeshPoints[0][0], MeshPoints[0][1], MeshPoints[0][2], # point in volume -> first box
                            0, # region tag (user-defined number)
                            1, # max tet volume in region
                            ]
    print "Building mesh from {} facets and {} points".format(len(MeshFacets), len(MeshPoints))
    try:
        mesh = build(mesh_info, options=Options(switches="pqT", epsilon=0.01), volume_constraints=True)
    except:
        pass
    try:
        mesh = build(mesh_info, options=Options(switches="pqT", epsilon=0.0001), volume_constraints=True)
    except:
        pass
    print "Created mesh with {} points, {} faces and {} elements.".format(len(mesh.points), len(mesh.faces), len(mesh.elements))

def create_mesh(big_r=1.0, small_r=0.5, num_points=10):
    dphi = 2 * np.pi / num_points

    # Compute the volume of a canonical tetrahedron
    # with edgelength radius2*dphi.
    a = small_r * dphi
    canonical_tet_volume = np.sqrt(2.0) / 12 * a ** 3

    radial_subdiv = int(2 * np.pi * big_r / a)

    rz = [
        (big_r + small_r * np.cos(i * dphi), 0.5 * small_r * np.sin(i * dphi))
        for i in range(num_points)
    ]

    geob = GeometryBuilder()
    geob.add_geometry(
        *generate_surface_of_revolution(
            rz, closure=EXT_CLOSED_IN_RZ, radial_subdiv=radial_subdiv
        )
    )
    mesh_info = MeshInfo()
    geob.set(mesh_info)
    meshpy_mesh = build(mesh_info, max_volume=canonical_tet_volume)

    return np.array(meshpy_mesh.points), np.array(meshpy_mesh.elements)

def make_cylinder_mesh(radius=0.5, height=1, radial_subdivisions=10,
        height_subdivisions=1, max_volume=None, periodic=False,
        boundary_tagger=(lambda fvi, el, fn, all_v: [])):
    from meshpy.tet import MeshInfo, build
    from meshpy.geometry import make_cylinder

    points, facets, facet_holestarts, facet_markers = \
            make_cylinder(radius, height, radial_subdivisions,
                    height_subdivisions)

    assert len(facets) == len(facet_markers)

    if periodic:
        return _make_z_periodic_mesh(
                points, facets, facet_holestarts, facet_markers,
                height=height,
                max_volume=max_volume,
                boundary_tagger=boundary_tagger)
    else:
        mesh_info = MeshInfo()
        mesh_info.set_points(points)
        mesh_info.set_facets_ex(facets, facet_holestarts, facet_markers)

        generated_mesh = build(mesh_info, max_volume=max_volume)

        from hedge.mesh import make_conformal_mesh_ext
        from hedge.mesh.element import Tetrahedron

        vertices = numpy.asarray(generated_mesh.points, dtype=float, order="C")
        return make_conformal_mesh_ext(
                vertices,
                [Tetrahedron(i, el_idx, vertices)
                    for i, el_idx in enumerate(generated_mesh.elements)],
                boundary_tagger)

def main():
    mesh_info = MeshInfo()

    mesh_info.set_points([
        (0, 0, 0),
        (2, 0, 0),
        (2, 2, 0),
        (0, 2, 0),
        (0, 0, 12),
        (2, 0, 12),
        (2, 2, 12),
        (0, 2, 12),
        ])

    mesh_info.set_facets([
        [0, 1, 2, 3],
        [4, 5, 6, 7],
        [0, 4, 5, 1],
        [1, 5, 6, 2],
        [2, 6, 7, 3],
        [3, 7, 4, 0],
        ])

    mesh_info.save_nodes("bar")
    mesh_info.save_poly("bar")

    mesh = build(mesh_info)

    mesh.save_nodes("barout")
    mesh.save_elements("barout")
    mesh.save_faces("barout")

    mesh.write_vtk("test.vtk")

def generateCubeMesh():
    mesh_info = MeshInfo()
    mesh_info.set_points([
        (0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0),
        (0, 0, 1), (1, 0, 1), (1, 1, 1), (0, 1, 1),
    ])
    mesh_info.set_facets([
        [0, 1, 2, 3],
        [4, 5, 6, 7],
        [0, 4, 5, 1],
        [1, 5, 6, 2],
        [2, 6, 7, 3],
        [3, 7, 4, 0],
    ])
    mesh = build(mesh_info)

    cellList = []
    vertexList = []
    mupifMesh = Mesh.UnstructuredMesh()

    print("Mesh Points:")
    for i, p in enumerate(mesh.points):
        print(i, p)
        vertexList.extend([Vertex.Vertex(i, i, p)])

    print("Point numbers in tetrahedra:")
    for i, t in enumerate(mesh.elements):
        print(i, t)
        cellList.extend([Cell.Tetrahedron_3d_lin(mupifMesh, i, i, t)])

    mupifMesh.setup(vertexList, cellList)

    return(mupifMesh)

def create_beam(vol):
    # build using MeshPy
    from meshpy.tet import MeshInfo,build

    mesh_info = MeshInfo()

    mesh_info.set_points([
        (0,0,0),
        (0,1,0),
        (0,1,1),
        (0,0,1),
        (5,0,0),
        (5,1,0),
        (5,1,1),
        (5,0,1),
        ])

    mesh_info.set_facets([
        [0,1,2,3],
        [4,5,6,7],
        [0,1,5,4],
        [1,2,6,5],
        [0,3,7,4],
        [3,2,6,7],
        ])


    mesh = build(mesh_info,max_volume=vol)

    return fmsh.MeshPyTet(mesh)

def main():
    from math import pi, cos, sin
    from meshpy.tet import MeshInfo, build
    from meshpy.geometry import generate_surface_of_revolution,\
            EXT_CLOSED_IN_RZ, GeometryBuilder

    big_r = 3
    little_r = 2.9

    points = 50
    dphi = 2*pi/points

    rz = [(big_r+little_r*cos(i*dphi), little_r*sin(i*dphi))
            for i in range(points)]

    geob = GeometryBuilder()
    geob.add_geometry(*generate_surface_of_revolution(rz,
            closure=EXT_CLOSED_IN_RZ, radial_subdiv=20))

    mesh_info = MeshInfo()
    geob.set(mesh_info)

    mesh_info.save_nodes("torus")
    mesh_info.save_poly("torus")
    mesh = build(mesh_info)
    mesh.write_vtk("torus.vtk")
    mesh.save_elements("torus_mesh")
    mesh.save_nodes("torus_mesh")

    mesh.write_neu(file("torus.neu", "w"),
            {1: ("pec", 0)})

def main():
    from math import pi, cos, sin
    from meshpy.tet import MeshInfo, build
    from meshpy.geometry import \
            generate_surface_of_revolution, EXT_OPEN, \
            GeometryBuilder

    r = 1
    l = 1

    rz = [(0,0), (r,0), (r,l), (0,l)]

    geob = GeometryBuilder()
    geob.add_geometry(*generate_surface_of_revolution(rz,
            radial_subdiv=20, ring_markers=[1,2,3]))

    mesh_info = MeshInfo()
    geob.set(mesh_info)

    mesh = build(mesh_info, max_volume=0.01)
    mesh.write_vtk("cylinder.vtk")
    mesh.write_neu(open("cylinder.neu", "w"), {
        1: ("minus_z", 1),
        2: ("outer", 2),
        3: ("plus_z", 3),
        })

def main():
    from math import pi, cos, sin
    from meshpy.tet import MeshInfo, build
    from meshpy.geometry import \
            generate_surface_of_revolution, EXT_OPEN, \
            GeometryBuilder

    r = 3

    points = 10
    dphi = pi/points

    def truncate(r):
        if abs(r) < 1e-10:
            return 0
        else:
            return r

    rz = [(truncate(r*sin(i*dphi)), r*cos(i*dphi)) for i in range(points+1)]

    geob = GeometryBuilder()
    geob.add_geometry(*generate_surface_of_revolution(rz,
            closure=EXT_OPEN, radial_subdiv=10))

    mesh_info = MeshInfo()
    geob.set(mesh_info)

    mesh = build(mesh_info)
    mesh.write_vtk("ball.vtk")

 def mesh(self):
     mesh_info = MeshInfo()
     mesh_info.set_points(transpose([self.flat_x,self.flat_y,self.flat_z]))
     mesh = build(mesh_info)
     print 'yo'
     for i, t in enumerate(mesh.elements):
         print i, t

def main():
    import numpy as np
    from math import pi, cos, sin
    from meshpy.tet import MeshInfo, build
    from meshpy.geometry import \
            GeometryBuilder, generate_surface_of_revolution, EXT_CLOSED_IN_RZ

    big_r = 3
    little_r = 1.5

    points = 50
    dphi = 2*pi/points

    rz = np.array([[big_r+little_r*cos(i*dphi), little_r*sin(i*dphi)]
            for i in range(points)])

    geo = GeometryBuilder()
    geo.add_geometry(
            *generate_surface_of_revolution(rz,
                closure=EXT_CLOSED_IN_RZ, radial_subdiv=20))

    mesh_info = MeshInfo()
    geo.set(mesh_info)

    mesh = build(mesh_info)

    def tet_face_vertices(vertices):
        return [(vertices[0], vertices[1], vertices[2]),
                (vertices[0], vertices[1], vertices[3]),
                (vertices[0], vertices[2], vertices[3]),
                (vertices[1], vertices[2], vertices[3]),
                ]

    face_map = {}
    for el_id, el in enumerate(mesh.elements):
        for fid, face_vertices in enumerate(tet_face_vertices(el)):
            face_map.setdefault(frozenset(face_vertices), []).append((el_id, fid))

    adjacency = {}
    for face_vertices, els_faces in face_map.items():
        if len(els_faces) == 2:
            (e1, f1), (e2, f2) = els_faces
            adjacency.setdefault(e1, []).append(e2)
            adjacency.setdefault(e2, []).append(e1)

    from pymetis import part_graph

    cuts, part_vert = part_graph(17, adjacency)

    try:
        import pyvtk
    except ImportError:
        print("Test succeeded, but could not import pyvtk to visualize result")
    else:
        vtkelements = pyvtk.VtkData(
            pyvtk.UnstructuredGrid(mesh.points, tetra=mesh.elements),
            "Mesh",
            pyvtk.CellData(pyvtk.Scalars(part_vert, name="partition")))
        vtkelements.tofile('split.vtk')

def test_tet_mesh(visualize=False):
    from math import pi, cos, sin
    from meshpy.tet import MeshInfo, build
    from meshpy.geometry import GeometryBuilder, generate_surface_of_revolution, EXT_CLOSED_IN_RZ

    pytest.importorskip("meshpy")

    big_r = 3
    little_r = 1.5

    points = 50
    dphi = 2*pi/points

    rz = np.array([[big_r+little_r*cos(i*dphi), little_r*sin(i*dphi)]
            for i in range(points)])

    geo = GeometryBuilder()
    geo.add_geometry(
            *generate_surface_of_revolution(rz,
                closure=EXT_CLOSED_IN_RZ, radial_subdiv=20))

    mesh_info = MeshInfo()
    geo.set(mesh_info)

    mesh = build(mesh_info)

    def tet_face_vertices(vertices):
        return [(vertices[0], vertices[1], vertices[2]),
                (vertices[0], vertices[1], vertices[3]),
                (vertices[0], vertices[2], vertices[3]),
                (vertices[1], vertices[2], vertices[3]),
                ]

    face_map = {}
    for el_id, el in enumerate(mesh.elements):
        for fid, face_vertices in enumerate(tet_face_vertices(el)):
            face_map.setdefault(frozenset(face_vertices), []).append((el_id, fid))

    adjacency = {}
    for face_vertices, els_faces in face_map.items():
        if len(els_faces) == 2:
            (e1, f1), (e2, f2) = els_faces
            adjacency.setdefault(e1, []).append(e2)
            adjacency.setdefault(e2, []).append(e1)

    import ipdb; ipdb.set_trace()
    cuts, part_vert = pymetis.part_graph(17, adjacency)

    if visualize:
        import pyvtk

        vtkelements = pyvtk.VtkData(
            pyvtk.UnstructuredGrid(mesh.points, tetra=mesh.elements),
            "Mesh",
            pyvtk.CellData(pyvtk.Scalars(part_vert, name="partition")))
        vtkelements.tofile('split.vtk')

def brep2lar(larBrep):
    V,FV = larBrep    
    mesh_info = MeshInfo()
    mesh_info.set_points(V)
    mesh_info.set_facets(FV)
    mesh = build(mesh_info)
    W = [v for h,v in enumerate(mesh.points)]
    CW = [tet for k,tet in enumerate(mesh.elements)]
    FW = sorted(set(AA(tuple)(CAT(AA(faces)(CW)))))
    EW = sorted(set(AA(tuple)(CAT(AA(edges)(FW)))))
    return W,CW,FW,EW

def test_tetgen_points():
    from meshpy.tet import MeshInfo, build, Options

    import numpy as np
    points = np.random.randn(10000, 3)

    mesh_info = MeshInfo()
    mesh_info.set_points(points)
    options = Options("")
    mesh = build(mesh_info, options=options)

    print(len(mesh.points))
    print(len(mesh.elements))

def test_torus():
    from math import pi, cos, sin
    from meshpy.tet import MeshInfo, build
    from meshpy.geometry import generate_surface_of_revolution, \
            EXT_CLOSED_IN_RZ, GeometryBuilder

    big_r = 3
    little_r = 2.9

    points = 50
    dphi = 2*pi/points

    rz = [(big_r+little_r*cos(i*dphi), little_r*sin(i*dphi))
            for i in range(points)]

    geob = GeometryBuilder()
    geob.add_geometry(*generate_surface_of_revolution(rz,
            closure=EXT_CLOSED_IN_RZ, radial_subdiv=20))

    mesh_info = MeshInfo()
    geob.set(mesh_info)

    build(mesh_info)

def _make_z_periodic_mesh(points, facets, facet_holestarts, facet_markers, height,
        max_volume, boundary_tagger):
    from meshpy.tet import MeshInfo, build
    from meshpy.geometry import Marker

    mesh_info = MeshInfo()
    mesh_info.set_points(points)
    mesh_info.set_facets_ex(facets, facet_holestarts, facet_markers)

    mesh_info.pbc_groups.resize(1)
    pbcg = mesh_info.pbc_groups[0]

    pbcg.facet_marker_1 = Marker.MINUS_Z
    pbcg.facet_marker_2 = Marker.PLUS_Z

    pbcg.set_transform(translation=[0, 0, height])

    def zper_boundary_tagger(fvi, el, fn, all_v):
        # we only ask about *boundaries*
        # we should not try to have the user tag
        # the (periodicity-induced) interior faces

        face_marker = fvi2fm[frozenset(fvi)]

        if face_marker == Marker.MINUS_Z:
            return ["minus_z"]
        if face_marker == Marker.PLUS_Z:
            return ["plus_z"]

        result = boundary_tagger(fvi, el, fn, all_v)
        if face_marker == Marker.SHELL:
            result.append("shell")
        return result

    generated_mesh = build(mesh_info, max_volume=max_volume)
    fvi2fm = generated_mesh.face_vertex_indices_to_face_marker

    from hedge.mesh import make_conformal_mesh_ext
    from hedge.mesh.element import Tetrahedron

    vertices = numpy.asarray(generated_mesh.points, dtype=float, order="C")
    return make_conformal_mesh_ext(
            vertices,
            [Tetrahedron(i, el_idx, vertices)
                for i, el_idx in enumerate(generated_mesh.elements)],
            zper_boundary_tagger,
            periodicity=[None, None, ("minus_z", "plus_z")])

def main():
    from meshpy.tet import MeshInfo, build
    from meshpy.geometry import generate_surface_of_revolution, GeometryBuilder

    simple_rz = [(0, 0), (1, 1), (1, 2), (0, 3)]

    geob = GeometryBuilder()
    geob.add_geometry(*generate_surface_of_revolution(simple_rz))

    mesh_info = MeshInfo()
    geob.set(mesh_info)

    # 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")

def main():
    from optparse import OptionParser

    description = "generate a vtk mesh from a superfish input file"
    parser = OptionParser(description=description)
    parser.add_option(
	    "--point-dist", dest="point_dist", default=0.3, type="float",
	    help="Spacing of intermediate points in circular arcs")
    parser.add_option(
	    "--inverse", dest="generate_inverse", action="store_true",
	    help="Whether to mesh the inverse of the gun (useless, but fun)")
    parser.add_option(
	    "--radial-subdiv", dest="radial_subdiv", action="store", type="int",
            default=16, help="How many radial subdivisions")
    parser.add_option(
	    "--tube-radius", dest="tube_radius", action="store", type="float",
            default=0.1, help="Radius of the interior tube")

    options, args = parser.parse_args()

    from superfish import parse_superfish_format
    rz = parse_superfish_format(args[0], options.point_dist)

    build_kwargs = {}

    if options.generate_inverse:
        from rzmesh import inverse_mesh_info_from_rz
        mesh_info = inverse_mesh_info_from_rz(rz, options.radial_subdiv)
    elif options.tube_radius:
        from rzmesh import make_mesh_info_with_inner_tube
        mesh_info = make_mesh_info_with_inner_tube(rz, 
                options.tube_radius, options.radial_subdiv)

        build_kwargs["volume_constraints"] = True
    else:
        from rzmesh import make_mesh_info
        mesh_info = make_mesh_info_with_inner_tube(rz, options.radial_subdiv)

    #mesh_info.save_poly("gun")
    #mesh_info.save_nodes("gun")

    from meshpy.tet import build
    mesh = build(mesh_info, verbose=True, **build_kwargs)
    print "%d elements" % len(mesh.elements)
    mesh.write_vtk("gun.vtk")

def main():
    from ply import parse_ply
    import sys
    data = parse_ply(sys.argv[1])

    from meshpy.geometry import GeometryBuilder

    builder = GeometryBuilder()
    builder.add_geometry(
            points=[pt[:3] for pt in data["vertex"].data],
            facets=[fd[0] for fd in data["face"].data])
    builder.wrap_in_box(1)

    from meshpy.tet import MeshInfo, build
    mi = MeshInfo()
    builder.set(mi)
    mi.set_holes([builder.center()])
    mesh = build(mi)
    print("%d elements" % len(mesh.elements))
    mesh.write_vtk("out.vtk")