How to use pygalmesh - 10 common examples

n1 = int(2 * numpy.pi * radius1 / self.edge_size)
            circ1 = [
                [
                    radius1 * numpy.cos((2 * numpy.pi * k) / n1),
                    radius1 * numpy.sin((2 * numpy.pi * k) / n1),
                    self.z1,
                ]
                for k in range(n1)
            ]
            circ1.append(circ1[0])
            return [circ0, circ1]

    edge_size = 0.12
    d = Hyperboloid(edge_size)

    mesh = pygalmesh.generate_mesh(d, cell_size=0.1, edge_size=edge_size, verbose=False)

    # TODO check the reference values
    tol = 1.0e-1
    assert abs(max(mesh.points[:, 0]) - 1.4) < tol
    assert abs(min(mesh.points[:, 0]) + 1.4) < tol
    assert abs(max(mesh.points[:, 1]) - 1.4) < tol
    assert abs(min(mesh.points[:, 1]) + 1.4) < tol
    assert abs(max(mesh.points[:, 2]) - 1.0) < tol
    assert abs(min(mesh.points[:, 2]) + 1.0) < tol

    vol = sum(helpers.compute_volumes(mesh.points, mesh.cells["tetra"]))
    assert abs(vol - 2 * numpy.pi * 47.0 / 60.0) < 0.16
    return

def test_torus():
    major_radius = 1.0
    minor_radius = 0.5
    s0 = pygalmesh.Torus(major_radius, minor_radius)
    mesh = pygalmesh.generate_mesh(s0, cell_size=0.1, verbose=False)

    tol = 1.0e-2
    radii_sum = major_radius + minor_radius
    assert abs(max(mesh.points[:, 0]) - radii_sum) < tol
    assert abs(min(mesh.points[:, 0]) + radii_sum) < tol
    assert abs(max(mesh.points[:, 1]) - radii_sum) < tol
    assert abs(min(mesh.points[:, 1]) + radii_sum) < tol
    assert abs(max(mesh.points[:, 2]) - minor_radius) < tol
    assert abs(min(mesh.points[:, 2]) + minor_radius) < tol

    vol = sum(helpers.compute_volumes(mesh.points, mesh.cells["tetra"]))
    ref_vol = (numpy.pi * minor_radius * minor_radius) * (2 * numpy.pi * major_radius)
    assert abs(vol - ref_vol) < 1.0e-1
    return

def generate_mesh():
    """Generates a fairly large mesh.
    """
    # import meshzoo
    # points, cells = meshzoo.rectangle(nx=300, ny=300)
    # return meshio.Mesh(points, {"triangle": cells})
    if os.path.isfile("cache.xdmf"):
        mesh = meshio.read("cache.xdmf")
    else:
        s = pygalmesh.Ball([0, 0, 0], 1.0)
        mesh = pygalmesh.generate_mesh(s, cell_size=2.0e-2, verbose=True)
        # mesh = pygalmesh.generate_mesh(s, cell_size=1.0e-1, verbose=True)
        mesh.cells = {"tetra": mesh.cells["tetra"]}
        mesh.point_data = []
        mesh.cell_data = {"tetra": {}}
        mesh.write("cache.xdmf")
    print(mesh)
    return mesh

c0 = pygalmesh.Cuboid([0, 0, -0.5], [3, 3, 0.5])
    c1 = pygalmesh.Cuboid([1, 1, -2], [2, 2, 2])
    u = pygalmesh.Intersection([c0, c1])

    # In CGAL, feature edges must not intersect, and that's a problem here: The
    # intersection edges of the cuboids share eight points with the edges of
    # the tall and skinny cuboid.
    # eps = 1.0e-2
    # extra_features = [
    #         [[1.0, 1.0 + eps, 0.5], [1.0, 2.0 - eps, 0.5]],
    #         [[1.0 + eps, 2.0, 0.5], [2.0 - eps, 2.0, 0.5]],
    #         [[2.0, 2.0 - eps, 0.5], [2.0, 1.0 + eps, 0.5]],
    #         [[2.0 - eps, 1.0, 0.5], [1.0 + eps, 1.0, 0.5]],
    #         ]

    mesh = pygalmesh.generate_mesh(u, cell_size=0.1, edge_size=0.1, verbose=False)

    # filter the vertices that belong to cells
    verts = mesh.points[numpy.unique(mesh.cells["tetra"])]

    tol = 1.0e-2
    assert abs(max(verts[:, 0]) - 2.0) < tol
    assert abs(min(verts[:, 0]) - 1.0) < tol
    assert abs(max(verts[:, 1]) - 2.0) < tol
    assert abs(min(verts[:, 1]) - 1.0) < tol
    assert abs(max(verts[:, 2]) - 0.5) < 0.05
    assert abs(min(verts[:, 2]) + 0.5) < 0.05

    vol = sum(helpers.compute_volumes(mesh.points, mesh.cells["tetra"]))
    assert abs(vol - 1.0) < 0.05

    return

radius = 1.0
    displacement = 0.5
    s0 = pygalmesh.Ball([displacement, 0, 0], radius)
    s1 = pygalmesh.Ball([-displacement, 0, 0], radius)
    u = pygalmesh.Difference(s0, s1)

    a = numpy.sqrt(radius ** 2 - displacement ** 2)
    edge_size = 0.15
    n = int(2 * numpy.pi * a / edge_size)
    circ = [
        [0.0, a * numpy.cos(i * 2 * numpy.pi / n), a * numpy.sin(i * 2 * numpy.pi / n)]
        for i in range(n)
    ]
    circ.append(circ[0])

    mesh = pygalmesh.generate_mesh(
        u,
        feature_edges=[circ],
        cell_size=0.15,
        edge_size=edge_size,
        facet_angle=25,
        facet_size=0.15,
        cell_radius_edge_ratio=2.0,
        verbose=False,
    )

    tol = 0.02
    assert abs(max(mesh.points[:, 0]) - (radius + displacement)) < tol
    assert abs(min(mesh.points[:, 0]) - 0.0) < tol
    assert abs(max(mesh.points[:, 1]) - radius) < tol
    assert abs(min(mesh.points[:, 1]) + radius) < tol
    assert abs(max(mesh.points[:, 2]) - radius) < tol

import os

import meshio
import matplotlib.pyplot as plt
import numpy
import pygalmesh

s = pygalmesh.Ball([0, 0, 0], 1.0)
mesh = pygalmesh.generate_mesh(s, cell_size=3.0e-2, verbose=True)
mesh.cells = {"tetra": mesh.cells["tetra"]}
mesh.point_data = {}
mesh.cell_data = {}

print("num points: {}".format(mesh.points.shape[0]))

formats = {
    "VTU (binary)": (lambda f, m: meshio.vtu.write(f, m, binary=True), ["out.vtu"]),
    "VTU (ASCII)": (lambda f, m: meshio.vtu.write(f, m, binary=False), ["out.vtu"]),
    "VTK (binary)": (lambda f, m: meshio.vtk.write(f, m, binary=True), ["out.vtk"]),
    "VTK (ASCII)": (lambda f, m: meshio.vtk.write(f, m, binary=False), ["out.vtk"]),
    "Gmsh 4.1 (binary)": (
        lambda f, m: meshio.gmsh.write(f, m, binary=True),
        ["out.msh"],
    ),
    "Gmsh 4.1 (ASCII)": (

def test_cone():
    base_radius = 1.0
    height = 2.0
    edge_size = 0.1
    s0 = pygalmesh.Cone(base_radius, height, edge_size)
    mesh = pygalmesh.generate_mesh(
        s0, cell_size=0.1, edge_size=edge_size, verbose=False
    )

    tol = 2.0e-1
    assert abs(max(mesh.points[:, 0]) - base_radius) < tol
    assert abs(min(mesh.points[:, 0]) + base_radius) < tol
    assert abs(max(mesh.points[:, 1]) - base_radius) < tol
    assert abs(min(mesh.points[:, 1]) + base_radius) < tol
    assert abs(max(mesh.points[:, 2]) - height) < tol
    assert abs(min(mesh.points[:, 2]) + 0.0) < tol

    vol = sum(helpers.compute_volumes(mesh.points, mesh.cells["tetra"]))
    ref_vol = numpy.pi * base_radius * base_radius / 3.0 * height
    assert abs(vol - ref_vol) < tol
    return

def test_ring_extrude():
    p = pygalmesh.Polygon2D([[0.5, -0.3], [1.5, -0.3], [1.0, 0.5]])
    edge_size = 0.1
    domain = pygalmesh.RingExtrude(p, edge_size)
    mesh = pygalmesh.generate_mesh(
        domain, cell_size=0.1, edge_size=edge_size, verbose=False
    )

    tol = 1.0e-2
    assert abs(max(mesh.points[:, 0]) - 1.5) < tol
    assert abs(min(mesh.points[:, 0]) + 1.5) < tol
    assert abs(max(mesh.points[:, 1]) - 1.5) < tol
    assert abs(min(mesh.points[:, 1]) + 1.5) < tol
    assert abs(max(mesh.points[:, 2]) - 0.5) < tol
    assert abs(min(mesh.points[:, 2]) + 0.3) < tol

    vol = sum(helpers.compute_volumes(mesh.points, mesh.cells["tetra"]))
    assert abs(vol - 2 * numpy.pi * 0.4) < 0.05
    return

def test_translation():
    s0 = pygalmesh.Translate(pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), [1.0, 0.0, 0.0])
    mesh = pygalmesh.generate_mesh(s0, cell_size=0.1, edge_size=0.1, verbose=False)

    tol = 1.0e-2
    assert abs(max(mesh.points[:, 0]) - 2.0) < tol
    assert abs(min(mesh.points[:, 0]) - 1.0) < tol
    assert abs(max(mesh.points[:, 1]) - 2.0) < tol
    assert abs(min(mesh.points[:, 1]) + 0.0) < tol
    assert abs(max(mesh.points[:, 2]) - 3.0) < tol
    assert abs(min(mesh.points[:, 2]) + 0.0) < tol
    vol = sum(helpers.compute_volumes(mesh.points, mesh.cells["tetra"]))
    assert abs(vol - 6.0) < tol
    return

def test_tetrahedron():
    s0 = pygalmesh.Tetrahedron(
        [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]
    )
    mesh = pygalmesh.generate_mesh(s0, cell_size=0.1, edge_size=0.1, verbose=False)

    tol = 1.0e-3
    assert abs(max(mesh.points[:, 0]) - 1.0) < tol
    assert abs(min(mesh.points[:, 0]) + 0.0) < tol
    assert abs(max(mesh.points[:, 1]) - 1.0) < tol
    assert abs(min(mesh.points[:, 1]) + 0.0) < tol
    assert abs(max(mesh.points[:, 2]) - 1.0) < tol
    assert abs(min(mesh.points[:, 2]) + 0.0) < tol

    vol = sum(helpers.compute_volumes(mesh.points, mesh.cells["tetra"]))
    assert abs(vol - 1.0 / 6.0) < tol
    return