def test_cli():
    input_mesh = helpers.tri_mesh
    infile = tempfile.NamedTemporaryFile().name
    meshio.write(infile, input_mesh, file_format="gmsh4-ascii")

    outfile = tempfile.NamedTemporaryFile().name

    meshio.cli.main(
        [
            infile,
            outfile,
            "--input-format",
            "gmsh-ascii",
            "--output-format",
            "vtk-binary",
        ]
    )

    mesh = meshio.read(outfile, file_format="vtk-binary")

    atol = 1.0e-15
    assert numpy.allclose(input_mesh.points, mesh.points, atol=atol, rtol=0.0)

    for cell_type, data in input_mesh.cells.items():
        assert numpy.allclose(data, mesh.cells[cell_type])
    return

def mesh(mesh_path, lc, H, L, b_dist, b_h, b_l, f_l, f_h, **NS_namespace):
    # Name mesh after local cell length
    mesh_path = mesh_path.format(lc)

    if not path.exists(mesh_path):
        # Initialize geometry
        geom = pygmsh.built_in.geometry.Geometry()

        # Surounding box
        b0 = geom.add_point([0, 0, 0], lcar=lc*2)
        b1 = geom.add_point([L, 0, 0], lcar=lc*2)
        b2 = geom.add_point([L, H, 0], lcar=lc*2)
        b3 = geom.add_point([0, H, 0], lcar=lc*2)

        # Inner geometry
        f0 = geom.add_point([b_dist, H / 2 - b_h / 2, 0], lcar=lc/3)
        f1 = geom.add_point([b_dist + b_l, H / 2 - b_h / 2, 0], lcar=lc/3)
        f2 = geom.add_point([b_dist + b_l, H / 2 - f_h / 2, 0], lcar=lc/3)
        f3 = geom.add_point([b_dist + b_l + f_l, H / 2 - f_h / 2, 0], lcar=lc/3)
        f4 = geom.add_point([b_dist + b_l + f_l, H / 2 + f_h / 2, 0], lcar=lc/3)
        f5 = geom.add_point([b_dist + b_l, H / 2 + f_h / 2, 0], lcar=lc/3)
        f6 = geom.add_point([b_dist + b_l, H / 2 + b_h / 2, 0], lcar=lc/3)
        f7 = geom.add_point([b_dist, H / 2 + b_h / 2, 0], lcar=lc/3)

        # Draw lines
        l0 = geom.add_line(b0, b1)
        l1 = geom.add_line(b1, b2)
        l2 = geom.add_line(b2, b3)
        l3 = geom.add_line(b3, b0)

        l4 = geom.add_line(f0, f1)
        l5 = geom.add_line(f1, f2)
        l6 = geom.add_line(f2, f3)
        l7 = geom.add_line(f3, f4)
        l8 = geom.add_line(f4, f5)
        l9 = geom.add_line(f5, f6)
        l10 = geom.add_line(f6, f7)
        l11 = geom.add_line(f7, f0)

        # Gather lines
        ll_outer = geom.add_line_loop(lines=[l0, l1, l2, l3])
        ll_inner = geom.add_line_loop(lines=[l4, l5, l6, l7, l8, l9, l10, l11])

        # Set surface
        ps = geom.add_plane_surface(ll_outer, [ll_inner])

        # Mesh surface
        points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom)

        # Write mesh
        meshio.write(mesh_path, meshio.Mesh(
                     points=points,
                     cells={"triangle": cells["triangle"]}))

    mesh = Mesh()
    with XDMFFile(MPI.comm_world, mesh_path) as infile:
        infile.read(mesh)

    return mesh

def mesh(mesh_path, mesh_function_path, lc, **NS_namespace):
    # Inspired from
    # https://gist.github.com/michalhabera/bbe8a17f788192e53fd758a67cbf3bed
    geom = pygmsh.built_in.geometry.Geometry()

    r1 = 1.7*2
    r2 = 3.8*2

    # Create geometry
    p0 = geom.add_point([ 0, 0, 0], lcar=lc)
    p1 = geom.add_point([ r1, 0, 0], lcar=lc)
    p2 = geom.add_point([ 0, -2*r2, 0], lcar=lc)
    p3 = geom.add_point([-r1, 0, 0], lcar=lc)

    p4 = geom.add_point([ 1.615, 0, 0], lcar=lc)
    p5 = geom.add_point([-0.085, 0, 0], lcar=lc)
    p6 = geom.add_point([-0.272, 0, 0], lcar=lc)
    p7 = geom.add_point([-1.632, 0, 0], lcar=lc)

    l0 = geom.add_ellipse_arc(p1, p0, p2, p2)
    l1 = geom.add_ellipse_arc(p2, p0, p3, p3)
    #l2 = geom.add_line(p3, p1)
    l3 = geom.add_line(p3, p7)
    l4 = geom.add_line(p7, p6)
    l5 = geom.add_line(p6, p5)
    l6 = geom.add_line(p5, p4)
    l7 = geom.add_line(p4, p1)

    #ll = geom.add_line_loop(lines=[l0, l1, l2])
    ll = geom.add_line_loop(lines=[l7, l0, l1, l3, l4, l5, l6])
    ps = geom.add_plane_surface(ll)

    # Tag line and surface
    geom.add_physical_line(lines=l4, label=1)
    geom.add_physical_line(lines=l6, label=2)
    geom.add_physical_line(lines=[l7, l0, l1, l3, l5], label=3)
    geom.add_physical_surface(surfaces=ps, label=4)

    # Mesh surface
    points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom)

    # Write, then read mesh and MeshFunction
    meshio.write(mesh_path, meshio.Mesh(
                 points=points,
                 cells={"triangle": cells["triangle"]}))

    meshio.write(mesh_function_path, meshio.Mesh(
                 points=points,
                 cells={"line": cells["line"]},
                 cell_data={"line": {"boundary": cell_data["line"]["gmsh:physical"]}}
                ))

    mesh = Mesh()
    with XDMFFile(MPI.comm_world, mesh_path) as infile:
        infile.read(mesh)

    return mesh

	def writeVTKFile(self,fn="",sub=False):
		
		"""Writes mesh into vtk file.
		
		Uses *meshIO* (https://github.com/nschloe/meshio), to convert the mesh saved
		in ``fnMesh`` to a .vtk file. 
		
		If ``sub==True``, will start a seperate subprocess and submit 
		``pyfrp_meshIO_script.py`` to it. This can be sometimes useful, since PyFRAP
		sometimes tends to crash otherwise.
		
		If no output path is given via ``fn``, will use same path as ``fnMesh``.
		
		.. note:: *meshIO* only gets imported inside this function, making PyFRAP running
		   even without the package installed. However, this feature will only run with
		   *meshIO*.
		
		Keyword Args:
			fn (str): Optional output path.
			sub (bool): Subprocess flag.
		
		Returns:
			str: Used output path.
		
		"""
		
		
		if not os.path.isfile(self.fnMesh):
			printWarning("Filepath to meshfile has not been specified yet. Cannot write VTK file.")
		
		if fn=="":
			fn=self.fnMesh.replace('.msh','.vtk')
		
		if sub:
			
			cmd = "python pyfrp_meshIO_script.py "+ self.fnMesh
			
			import shlex
			import subprocess
			args = shlex.split(cmd)
			
			p = subprocess.Popen(args)
			p.wait()
		else:
			#MeshIO
			import meshio
			
			points, cells, point_data, cell_data, field_data = meshio.read(self.fnMesh)
			meshio.write(fn,points,cells,point_data=point_data,cell_data=cell_data,field_data=field_data)
			
		return fn

    def gen_mesh(self, filename, n):
        import pygmsh as pg
        geom = pg.Geometry()

        loops = []
        for be in self.blade_elements:
            car = 0.5/1000
            line_l, line_u = be.get_foil_points(n, self.get_scimitar_offset(be.r))
            loop_points = np.concatenate((line_l, line_u[::-1]), axis=0)
            g_pts = []
            for p in loop_points[0:-2]:
              g_pts.append(geom.add_point(p,car))
              
            l_foil = geom.add_bspline(g_pts)
            
            loops.append(l_foil)

        #print geom.get_code()
        geom.add_ruled_surface([loops[0], loops[1]])
        
        #l in range(len(loops) - 1):
            #geom.add_surface(loops[l], loops[l+1])
        #poly = geom.add_polygon([
            #[0.0,   0.5, 0.0],
            #[-0.1,  0.1, 0.0],
            #[-0.5,  0.0, 0.0],
            #[-0.1, -0.1, 0.0],
            #[0.0,  -0.5, 0.0],
            #[0.1,  -0.1, 0.0],
            #[0.5,   0.0, 0.0],
            #[0.1,   0.1, 0.0]
            #],
            #lcar=0.05
            #)
        #axis = [0, 0, 1]

        #geom.extrude(
            #'Surface{%s}' % poly,
            #translation_axis=axis,
            #rotation_axis=axis,
            #point_on_axis=[0, 0, 0],
            #angle=2.0 / 6.0 * np.pi
            #)

        points, cells = pg.generate_mesh(geom)

        import meshio
        meshio.write(filename, points, cells)

def write_file(name, obj, **cell_data):

    for ext in ['msh','vtk']:
        t1 = time.time()
        fname = "test_geometry_%s.%s"%(name, ext)

        print 'min/max elements:', numpy.min(obj.grid_elements), numpy.max(obj.grid_elements)
        print "points: ",obj.grid_points.shape

        cd = dict(cell_data, domains=obj.grid_domains)
        print "cell data:", obj.grid_elements.shape
        for k,v in cd.items():
            print "\t",k,v.shape
        meshio.write(fname, obj.grid_points,
                     dict(triangle=obj.grid_elements),
                     cell_data = cd)
        t2 = time.time()
        print 'wrote %s in %.1f seconds' % (fname, t2-t1)

def _write_read(filename, mesh):
    '''Write and read a file, and make sure the data is the same as before.
    '''
    try:
        input_point_data = mesh['point_data']
    except KeyError:
        input_point_data = {}

    try:
        input_cell_data = mesh['cell_data']
    except KeyError:
        input_cell_data = {}

    meshio.write(
        filename,
        mesh['points'], mesh['cells'],
        point_data=input_point_data,
        cell_data=input_cell_data
        )
    points, cells, point_data, cell_data, _ = meshio.read(filename)

    # Numpy's array_equal is too strict here, cf.
    # <https://mail.scipy.org/pipermail/numpy-discussion/2015-December/074410.html>.
    # Use allclose.

    # We cannot compare the exact rows here since the order of the points might
    # have changes. Just compare the sums
    assert numpy.allclose(mesh['points'], points)

    for key, data in mesh['cells'].items():
        assert numpy.allclose(data, cells[key])
    for key, data in input_point_data.items():
        assert numpy.allclose(data, point_data[key])
    for key, data in input_cell_data.items():
        assert numpy.allclose(data, cell_data[key])

    os.remove(filename)

    return

def _write_read(filename, mesh):
    '''Write and read a file, and make sure the data is the same as before.
    '''
    meshio.write(
        filename,
        mesh['points'], mesh['cells'],
        point_data=mesh['point_data'],
        cell_data=mesh['cell_data']
        )
    points, cells, point_data, cell_data, _ = meshio.read(filename)

    # We cannot compare the exact rows here since the order of the points might
    # have changes. Just compare the sums
    assert numpy.allclose(mesh['points'], points)

    for key, data in mesh['cells'].items():
        assert numpy.array_equal(data, cells[key])
    for key, data in mesh['point_data'].items():
        assert numpy.array_equal(data, point_data[key])
    for key, data in mesh['cell_data'].items():
        assert numpy.array_equal(data, cell_data[key])
    return

def read_write():
    X, cells = generate_mesh()

    formats = [
        "ansys-ascii",
        "ansys-binary",
        "exodus",
        "dolfin-xml",
        "gmsh-ascii",
        "gmsh-binary",
        "med",
        "medit",
        "permas",
        "moab",
        "off",
        "stl-ascii",
        "stl-binary",
        "vtk-ascii",
        "vtk-binary",
        "vtu-ascii",
        "vtu-binary",
        "xdmf",
    ]

    filename = "foo"
    print()
    print("format        write (s)    read(s)")
    print()
    for fmt in formats:
        t = time.time()
        meshio.write(filename, X, cells, file_format=fmt)
        elapsed_write = time.time() - t

        t = time.time()
        meshio.read(filename, file_format=fmt)
        elapsed_read = time.time() - t
        print("{0: <12}  {1:e} {2:e}".format(fmt, elapsed_write, elapsed_read))

    return

def _write_file(name, p, e, pd=None, cd=None, fd=None):
    pd = pd or dict()
    cd = cd or dict()
    fd = fd or dict()

    print 'points:',p.shape
    print 'elements:',len(e)
    for k,v in e.items():
        print '\t%s:%s' % (k, v.shape)
    print 'point data:',len(pd)
    for k,v in pd.items():
        print '\t%s:%s' % (k,v.shape)
    print 'cell data:',len(cd)
    for k,v in cd.items():
        print '\t%s:%s' % (k,v.shape)

    for ext in ['msh','vtk']:
        t1 = time.time()
        fname = "test_geometry_%s.%s"%(name, ext)
        meshio.write(fname,p,e,
                     point_data=pd,cell_data=cd,field_data=fd)
        t2 = time.time()
        print 'wrote %s in %.1f seconds' % (fname, t2-t1)

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


if __name__ == "__main__":
    import meshio

    points, cells = create_mesh()
    meshio.write("torus.e", points, {"tetra": cells})

        (0.5 * lx, 0.0, 0.5 * lz),
        (0.0, 0.0, 0.5 * lz),
        (0.0, 0.0, 0.0),
    ]
    facets = [
        [0, 1, 2, 3],
        [4, 7, 12, 11, 5, 6],
        [0, 1, 9, 8, 7, 4],
        [1, 2, 5, 11, 10, 9],
        [2, 5, 6, 3],
        [3, 6, 4, 0],
        [8, 13, 12, 7],
        [8, 9, 10, 13],
        [10, 11, 12, 13],
    ]
    # create the mesh
    # Set the geometry and build the mesh.
    info = meshpy.tet.MeshInfo()
    info.set_points(points)
    info.set_facets(facets)
    meshpy_mesh = meshpy.tet.build(info, max_volume=maxvol)

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


if __name__ == "__main__":
    import meshio

    points, cells = create_mesh()
    meshio.write("lshape3d.e", points, {"tetra": cells})

# -*- coding: utf-8 -*-

import pygmsh


def test(lcar=1.0):
    geom = pygmsh.built_in.Geometry()
    poly = geom.add_polygon(
        [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0]], lcar
    )

    geom.set_transfinite_surface(poly.surface, size=[11, 9])

    mesh = pygmsh.generate_mesh(geom, geo_filename="transfinite.geo")
    assert len(mesh.cells["triangle"]) == 10 * 8 * 2
    return mesh


if __name__ == "__main__":
    import meshio

    meshio.write("transfinite.vtu", test())

# -*- coding: utf-8 -*-
import pygmsh

from helpers import compute_volume


def test():
    geom = pygmsh.built_in.Geometry()

    geom.add_circle(
        [0.0, 0.0, 0.0],
        1.0,
        lcar=0.1,
        num_sections=4,
        # If compound==False, the section borders have to be points of the
        # discretization. If using a compound circle, they don't; gmsh can
        # choose by itself where to point the circle points.
        compound=True,
    )

    ref = 3.1363871677682247
    mesh = pygmsh.generate_mesh(geom, prune_z_0=True)
    assert abs(compute_volume(mesh) - ref) < 1.0e-2 * ref
    return mesh


if __name__ == "__main__":
    import meshio

    meshio.write("circle.vtk", test())

# -*- coding: utf-8 -*-
"""
Creates a mesh for a square with a round hole.
"""
import pygmsh

from helpers import compute_volume


def test():
    geom = pygmsh.built_in.Geometry()

    circle = geom.add_circle(
        x0=[0.5, 0.5, 0.0], radius=0.25, lcar=0.1, num_sections=4, make_surface=False
    )

    geom.add_rectangle(0.0, 1.0, 0.0, 1.0, 0.0, lcar=0.1, holes=[circle.line_loop])

    ref = 0.8086582838174551
    mesh = pygmsh.generate_mesh(geom, geo_filename="h.geo")
    assert abs(compute_volume(mesh) - ref) < 1.0e-2 * ref
    return mesh


if __name__ == "__main__":
    import meshio

    meshio.write("rectangle_with_hole.vtu", test())

    geom = pygmsh.built_in.Geometry()

    # Draw a square
    poly = geom.add_polygon(
        [[+0.5, +0.0, 0.0], [+0.0, +0.5, 0.0], [-0.5, +0.0, 0.0], [+0.0, -0.5, 0.0]],
        lcar=lcar,
    )

    axis = [0, 0, 1.0]

    geom.extrude(
        poly,
        translation_axis=axis,
        rotation_axis=axis,
        point_on_axis=[0.0, 0.0, 0.0],
        angle=0.5 * pi,
        num_layers=5,
        recombine=True,
    )

    ref = 3.98156496566
    mesh = pygmsh.generate_mesh(geom)
    assert abs(compute_volume(mesh) - ref) < 1.0e-2 * ref
    return mesh


if __name__ == "__main__":
    import meshio

    meshio.write("rotated_layers.vtu", test())

        )
    )

    for alpha, a1, z in zip(Alpha, A1, Z_pos):
        # Rotate torus to the y-z-plane.
        R1 = pygmsh.rotation_matrix([0.0, 1.0, 0.0], 0.5 * np.pi)
        R2 = pygmsh.rotation_matrix([0.0, 0.0, 1.0], alpha)
        x0 = np.array([a1, 0.0, 0.0])
        x1 = np.array([0.0, 0.0, z])
        # First rotate to y-z-plane, then move out to a1, rotate by angle
        # alpha, move up by z.
        #
        #    xnew = R2*(R1*x+x0) + x1
        #
        geom.add_torus(
            irad=irad, orad=orad, lcar=0.1, R=np.dot(R2, R1), x0=np.dot(R2, x0) + x1
        )

    geom.add_box(-1.0, 1.0, -1.0, 1.0, -1.0, 1.0, lcar=0.3)

    ref = len(A1) * 2 * np.pi ** 2 * orad * irad ** 2 + 2.0 ** 3
    mesh = pygmsh.generate_mesh(geom)
    assert np.isclose(compute_volume(mesh), ref, rtol=2e-2)
    return mesh


if __name__ == "__main__":
    import meshio

    meshio.write("torus_crowd.vtu", test())

        R2 = pg.rotation_matrix([0.0, 0.0, 1.0], alpha)
        x0 = np.array([a1, 0.0, 0.0])
        x1 = np.array([0.0, 0.0, z])
        # First rotate to y-z-plane, then move out to a1, rotate by angle
        # alpha, move up by z.
        #
        #    xnew = R2*(R1*x+x0) + x1
        #
        geom.add_torus(
                irad=irad,
                orad=orad,
                lcar=0.1,
                R=np.dot(R2, R1),
                x0=np.dot(R2, x0) + x1
                )

    geom.add_box(
            -1.0, 1.0,
            -1.0, 1.0,
            -1.0, 1.0,
            lcar=0.3
            )

    return geom


if __name__ == '__main__':
    import meshio
    points, cells = pg.generate_mesh(generate())
    meshio.write('torus_crowd.vtu', points, cells)

# -*- coding: utf-8 -*-

import pygmsh


def test(lcar=0.5):
    geom = pygmsh.built_in.Geometry()
    poly = geom.add_polygon(
        [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0]], lcar
    )

    top, volume, lat = geom.extrude(poly.surface, [0, 0, 2])

    geom.add_physical(poly.surface, label="bottom")
    geom.add_physical(top, label="top")
    geom.add_physical(volume, label="volume")
    geom.add_physical(lat, label="lat")
    geom.add_physical(poly.lines[0], label="line")

    mesh = pygmsh.generate_mesh(geom, geo_filename="physical.geo")
    return mesh


if __name__ == "__main__":
    import meshio

    meshio.write("physical.vtu", test())

    # Create polygon for airfoil
    char_length = 1.0e-1
    airfoil = geom.add_polygon(airfoil_coordinates, char_length, make_surface=False)

    # Create surface for numerical domain with an airfoil-shaped hole
    left_dist = 1.0
    right_dist = 3.0
    top_dist = 1.0
    bottom_dist = 1.0
    xmin = airfoil_coordinates[:, 0].min() - left_dist * coord
    xmax = airfoil_coordinates[:, 0].max() + right_dist * coord
    ymin = airfoil_coordinates[:, 1].min() - bottom_dist * coord
    ymax = airfoil_coordinates[:, 1].max() + top_dist * coord
    domainCoordinates = numpy.array(
        [[xmin, ymin, 0.0], [xmax, ymin, 0.0], [xmax, ymax, 0.0], [xmin, ymax, 0.0]]
    )
    polygon = geom.add_polygon(domainCoordinates, char_length, holes=[airfoil])
    geom.add_raw_code("Recombine Surface {{{}}};".format(polygon.surface.id))

    ref = 10.525891646546
    mesh = pygmsh.generate_mesh(geom, remove_faces=True)
    assert abs(compute_volume(mesh) - ref) < 1.0e-2 * ref
    return mesh


if __name__ == "__main__":
    import meshio

    meshio.write("airfoil.vtu", test())