def test():
    geom = pygmsh.built_in.Geometry()
    geom.add_circle([0, 0, 0], 1, 0.1, make_surface=False)
    mesh = pygmsh.generate_mesh(geom)
    ref = 2 * np.pi
    assert np.abs(compute_volume(mesh) - ref) < 1e-2 * ref
    return

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

    # Draw a cross with a circular hole
    circ = geom.add_circle([0.0, 0.0, 0.0], 0.1, lcar=lcar)
    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=lcar,
        holes=[circ],
    )

    axis = [0, 0, 1.0]

    geom.extrude(
        poly,
        translation_axis=axis,
        rotation_axis=axis,
        point_on_axis=[0, 0, 0],
        angle=2.0 / 6.0 * np.pi,
    )

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

def test():
    geom = pygmsh.opencascade.Geometry(
        characteristic_length_min=2.0, characteristic_length_max=2.0
    )

    rect1 = geom.add_rectangle([10.0, 0.0, 0.0], 20.0, 40.0, corner_radius=5.0)
    rect2 = geom.add_rectangle([0.0, 10.0, 0.0], 40.0, 20.0, corner_radius=5.0)
    disk1 = geom.add_disk([14.5, 35.0, 0.0], 1.85)
    disk2 = geom.add_disk([25.5, 5.0, 0.0], 1.85)

    rect3 = geom.add_rectangle([10.0, 30.0, 0.0], 10.0, 1.0)
    rect4 = geom.add_rectangle([20.0, 9.0, 0.0], 10.0, 1.0)

    r1 = geom.add_rectangle([9.0, 0.0, 0.0], 21.0, 20.5, corner_radius=8.0)
    r2 = geom.add_rectangle([10.0, 00.0, 0.0], 20.0, 19.5, corner_radius=7.0)
    diff1 = geom.boolean_difference([r1], [r2])
    r22 = geom.add_rectangle([9.0, 10.0, 0.0], 11.0, 11.0)
    inter1 = geom.boolean_intersection([diff1, r22])

    r3 = geom.add_rectangle([10.0, 19.5, 0.0], 21.0, 21.0, corner_radius=8.0)
    r4 = geom.add_rectangle([10.0, 20.5, 0.0], 20.0, 20.0, corner_radius=7.0)
    diff2 = geom.boolean_difference([r3], [r4])
    r33 = geom.add_rectangle([20.0, 19.0, 0.0], 11.0, 11.0)
    inter2 = geom.boolean_intersection([diff2, r33])

    geom.boolean_difference(
        [rect1, rect2], [disk1, disk2, rect3, rect4, inter1, inter2]
    )

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

def test(lcar=1.0):
    lbw = [2, 3, 5]
    geom = pygmsh.built_in.Geometry()
    points = [geom.add_point([x, 0.0, 0.0], lcar) for x in [0.0, lbw[0]]]
    line = geom.add_line(*points)

    _, rectangle, _ = geom.extrude(
        line, translation_axis=[0.0, lbw[1], 0.0], num_layers=lbw[1], recombine=True
    )
    geom.extrude(
        rectangle,
        translation_axis=[0.0, 0.0, lbw[2]],
        num_layers=lbw[2],
        recombine=True,
    )

    # compute_volume only supports 3D for tetras, but does return
    # surface area for quads

    ref = sum(l * w for l, w in permutations(lbw, 2))  # surface area
    mesh = pygmsh.generate_mesh(geom, prune_vertices=False)
    # TODO compute hex volumes
    assert (
        abs(
            compute_volume(meshio.Mesh(mesh.points, {"quad": mesh.cells["quad"]})) - ref
        )
        < 1.0e-2 * ref
    )
    return mesh

def test_square_circle_slice():
    """Test planar suface square with circular hole.

    Also test for surface area of fragments.
    """
    geo_object = built_in_opencascade_geos_fragments()

    # Gmsh 4 default format MSH4 doesn't have geometrical entities.
    mesh = pygmsh.generate_mesh(geo_object, extra_gmsh_arguments=["-format", "msh2"])
    ref = 1
    val = compute_volume(mesh)
    assert np.abs(val - ref) < 1e-3 * ref

    outer_mask = np.where(mesh.cell_data["triangle"]["gmsh:geometrical"] == 13)[0]
    outer_cells = {}
    outer_cells["triangle"] = mesh.cells["triangle"][outer_mask]

    inner_mask = np.where(mesh.cell_data["triangle"]["gmsh:geometrical"] == 12)[0]
    inner_cells = {}
    inner_cells["triangle"] = mesh.cells["triangle"][inner_mask]

    ref = 1 - 0.1 ** 2 * np.pi
    value = compute_volume(meshio.Mesh(mesh.points, outer_cells))
    assert np.abs(value - ref) < 1e-2 * ref
    return

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

    X0 = np.array(
        [[+0.0, +0.0, 0.0], [+0.5, +0.3, 0.1], [-0.5, +0.3, 0.1], [+0.5, -0.3, 0.1]]
    )

    R = np.array([0.1, 0.2, 0.1, 0.14])

    holes = [
        geom.add_ball(x0, r, with_volume=False, lcar=0.2 * r).surface_loop
        for x0, r in zip(X0, R)
    ]

    # geom.add_box(
    #         -1, 1,
    #         -1, 1,
    #         -1, 1,
    #         lcar=0.2,
    #         holes=holes
    #         )

    geom.add_ball([0, 0, 0], 1.0, lcar=0.2, holes=holes)

    # geom.add_physical_volume(ball, label="cheese")

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

def test_fragments_diff_union():
    """Test planar surface with holes.

    Construct it with boolean operations and verify that it is the same.
    """
    # construct surface using boolean
    geo_object = pygmsh.opencascade.Geometry(0.04, 0.04)
    geo_object, square = square_loop(geo_object)
    geo_object, line_loop = circle_loop(geo_object)
    surf1 = geo_object.add_plane_surface(square)
    surf2 = geo_object.add_plane_surface(line_loop)

    geo_object.add_physical([surf1], label=1)
    geo_object.add_physical([surf2], label=2)
    surf_diff = geo_object.boolean_difference([surf1], [surf2], delete_other=False)
    geo_object.boolean_union([surf_diff, surf2])
    mesh = pygmsh.generate_mesh(geo_object)
    assert np.abs((compute_volume(mesh) - 1) / 1) < 1e-3
    surf = 1 - 0.1 ** 2 * np.pi
    outer_mask = np.where(mesh.cell_data["triangle"]["gmsh:physical"] == 1)[0]
    outer_cells = {}
    outer_cells["triangle"] = mesh.cells["triangle"][outer_mask]

    inner_mask = np.where(mesh.cell_data["triangle"]["gmsh:physical"] == 2)[0]
    inner_cells = {}
    inner_cells["triangle"] = mesh.cells["triangle"][inner_mask]

    value = compute_volume(meshio.Mesh(mesh.points, outer_cells))
    assert np.abs((value - surf)) < 1e-2 * surf
    return

def test(lcar=0.05):
    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

def test():
    # Characteristic length
    lcar = 1e-1

    # Coordinates of lower-left and upper-right vertices of a square domain
    xmin = 0.0
    xmax = 5.0
    ymin = 0.0
    ymax = 5.0

    # Vertices of a square hole
    squareHoleCoordinates = np.array([[1, 1, 0], [4, 1, 0], [4, 4, 0], [1, 4, 0]])

    # Create geometric object
    geom = pygmsh.built_in.Geometry()

    # Create square hole
    squareHole = geom.add_polygon(squareHoleCoordinates, lcar, make_surface=False)

    # Create square domain with square hole
    geom.add_rectangle(xmin, xmax, ymin, ymax, 0.0, lcar, holes=[squareHole.line_loop])

    mesh = pygmsh.generate_mesh(geom, extra_gmsh_arguments=["-order", "2"])
    # TODO support for volumes of triangle6
    # ref = 16.0
    # from helpers import compute_volume
    # assert abs(compute_volume(points, cells) - ref) < 1.0e-2 * ref
    return mesh

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

    # Draw a cross with a circular hole
    circ = geom.add_circle([0.0, 0.0, 0.0], 0.1, lcar=lcar, make_surface=False)
    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=lcar,
        holes=[circ],
    )

    axis = [0, 0, 1.0]

    geom.extrude(poly, translation_axis=axis, num_layers=1)

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

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 test():
    geom = pygmsh.built_in.Geometry()
    geom.add_ellipsoid([0.0, 0.0, 0.0], [1.0, 0.5, 0.75], 0.05)

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

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

    geom.add_torus([0.0, 0.0, 0.0], 1.0, 0.3, 1.25 * pi, char_length=0.1)

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

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

    geom.add_cylinder([0.0, 0.0, 0.0], [0.0, 0.0, 1.0], 0.5, 0.25 * pi, char_length=0.1)

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

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

    geom.add_rectangle(0.0, 1.0, 0.0, 1.0, 0.0, 0.1)

    ref = 1.0
    mesh = pygmsh.generate_mesh(geom, mesh_file_type="vtk")
    assert abs(compute_volume(mesh) - ref) < 1.0e-2 * ref
    return mesh

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

    geom.add_wedge([0.0, 0.0, 0.0], [1.0, 1.0, 1.0], top_extent=0.4, char_length=0.1)

    ref = 0.7
    mesh = pygmsh.generate_mesh(geom)
    assert abs(compute_volume(mesh) - ref) < 1.0e-2 * ref
    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 test():
    geom = pygmsh.opencascade.Geometry()

    geom.add_box([0.0, 0.0, 0.0], [1, 2, 3], char_length=0.1)

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

def test_square_circle_hole():
    """Test planar surface with holes.

    Construct it with boolean operations and verify that it is the same.
    """
    for geo_object in built_in_opencascade_geos():
        mesh = pygmsh.generate_mesh(geo_object)
        surf = 1 - 0.1 ** 2 * np.pi
        assert np.abs((compute_volume(mesh) - surf) / surf) < 1e-3
    return

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

    geom.add_cone(
        [0.0, 0.0, 0.0], [0.0, 0.0, 1.0], 1.0, 0.3, 1.25 * pi, char_length=0.1
    )

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