Helper Functions and Utilities

def write(filename: str):
    """
    Write current Gmsh model to file.
    
    Supports multiple file formats based on extension:
    - .msh: Gmsh mesh format (versions 2.2, 4.0, 4.1)
    - .vtk: VTK format for ParaView/VisIt visualization
    - .vtu: VTK unstructured grid format
    - .med: MED format for Code_Aster/Salome
    - .stl: STL format for 3D printing
    - .ply: PLY format for 3D graphics
    - .obj: Wavefront OBJ format
    - .off: Object File Format
    
    Parameters:
    - filename: str, output filename with extension determining format
    
    Returns:
    None (writes file to disk)
    """

def rotation_matrix(u, theta):
    """
    Create rotation matrix for rotation around vector u by angle theta.
    
    Based on Rodrigues' rotation formula:
    https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
    
    Parameters:
    - u: array-like, rotation axis vector (must be unit vector)
    - theta: float, rotation angle in radians
    
    Returns:
    numpy.ndarray: 3x3 rotation matrix
    
    Raises:
    AssertionError: if u is not a unit vector (|u| != 1)
    """

def orient_lines(lines):
    """
    Reorder and reorient unordered lines to form a closed polygon.
    
    Takes a sequence of unordered and potentially mis-oriented line segments
    that define a closed polygon and returns them in proper order with 
    correct orientations to form a valid curve loop.
    
    Parameters:
    - lines: sequence of Line objects defining a closed polygon
    
    Returns:
    list: Reordered and reoriented Line objects forming closed loop
    
    Algorithm:
    1. Analyzes vertex connectivity of all line segments
    2. Reorders lines to form continuous path
    3. Reverses line orientation where needed for consistency
    4. Ensures resulting loop is closed and properly oriented
    """

def optimize(mesh, method="", verbose=False):
    """
    Optimize mesh quality using Gmsh optimization methods.
    
    Applies mesh optimization algorithms to improve element quality metrics
    such as aspect ratio, skewness, and Jacobian determinant. Optimization
    can significantly improve numerical accuracy in finite element analysis.
    
    Parameters:
    - mesh: meshio.Mesh, input mesh object to optimize
    - method: str, optimization method name (empty string uses default)
      Available methods include:
      - "": Default optimization (Gmsh automatic selection)
      - "Netgen": Netgen-based optimization
      - "HighOrder": High-order node positioning optimization
      - "HighOrderElastic": Elastic high-order optimization
      - "HighOrderFastCurving": Fast curving for high-order elements
      - "Laplace2D": Laplacian smoothing for 2D meshes
      - "Relocate2D": Node relocation for 2D meshes
    - verbose: bool, enable verbose output during optimization
    
    Returns:
    meshio.Mesh: Optimized mesh object with improved element quality
    
    Note:
    Optimization modifies node positions and may add/remove elements
    to improve overall mesh quality. Physical groups and mesh regions
    are preserved during the optimization process.
    """

def extract_to_meshio():
    """
    Extract current Gmsh mesh to meshio format.
    
    Converts the active Gmsh mesh model into a meshio.Mesh object,
    including all mesh entities (points, lines, surfaces, volumes),
    physical groups, and material properties.
    
    Returns:
    meshio.Mesh: Mesh object containing:
    - points: numpy.ndarray, node coordinates (N×3)
    - cells: list of CellBlock objects with connectivity data
    - cell_sets: dict, physical groups mapped to cell indices
    - point_data: dict, node-based data fields (if any)
    - cell_data: dict, element-based data fields (if any)
    
    Supported element types:
    - 1D: Line elements (line, line3)
    - 2D: Triangle, quadrilateral elements (triangle, triangle6, quad, quad9)
    - 3D: Tetrahedron, hexahedron, prism, pyramid elements
    
    Physical groups are converted to cell_sets with original names,
    enabling material assignment and boundary condition specification
    in external analysis codes.
    """

import pygmsh

# Create geometry and generate initial mesh
with pygmsh.occ.Geometry() as geom:
    # Create complex geometry
    sphere = geom.add_ball([0, 0, 0], 1.0, mesh_size=0.2)
    box = geom.add_box([-0.5, -0.5, -0.5], [1, 1, 1], mesh_size=0.15)
    complex_shape = geom.boolean_union([sphere, box])
    
    # Generate initial mesh
    initial_mesh = geom.generate_mesh(dim=3)

# Optimize mesh quality
optimized_mesh = pygmsh.optimize(initial_mesh, method="HighOrder", verbose=True)

# Write optimized mesh to file
pygmsh.write("optimized_mesh.msh")

# Compare mesh statistics
print(f"Initial elements: {len(initial_mesh.cells[0].data)}")
print(f"Optimized elements: {len(optimized_mesh.cells[0].data)}")

import pygmsh
import numpy as np

with pygmsh.geo.Geometry() as geom:
    # Create unordered line segments for a polygon
    points = [
        geom.add_point([0, 0, 0]),
        geom.add_point([1, 0, 0]),
        geom.add_point([1, 1, 0]),
        geom.add_point([0.5, 1.5, 0]),
        geom.add_point([0, 1, 0])
    ]
    
    # Create lines in random order (simulating imported data)
    unordered_lines = [
        geom.add_line(points[2], points[3]),  # line 2-3
        geom.add_line(points[0], points[1]),  # line 0-1
        geom.add_line(points[4], points[0]),  # line 4-0
        geom.add_line(points[1], points[2]),  # line 1-2
        geom.add_line(points[3], points[4])   # line 3-4
    ]
    
    # Orient lines to form proper closed loop
    oriented_lines = pygmsh.orient_lines(unordered_lines)
    
    # Create surface from properly oriented lines
    curve_loop = geom.add_curve_loop(oriented_lines)
    surface = geom.add_plane_surface(curve_loop)
    
    mesh = geom.generate_mesh(dim=2)

import pygmsh
import numpy as np

with pygmsh.occ.Geometry() as geom:
    # Create base geometry
    cylinder = geom.add_cylinder([0, 0, 0], [0, 0, 1], 0.5)
    
    # Create rotation matrix for 45-degree rotation around [1,1,1] axis
    axis = np.array([1, 1, 1]) / np.sqrt(3)  # Normalize to unit vector
    angle = np.pi / 4  # 45 degrees
    R = pygmsh.rotation_matrix(axis, angle)
    
    # Apply rotation to create multiple rotated copies
    rotated_cylinders = [cylinder]
    for i in range(1, 8):
        # Accumulate rotations
        total_angle = i * angle
        R_total = pygmsh.rotation_matrix(axis, total_angle)
        
        # Create rotated copy (conceptual - actual rotation done via Gmsh)
        rotated_cyl = geom.copy(cylinder)
        geom.rotate(rotated_cyl, [0, 0, 0], total_angle, axis)
        rotated_cylinders.append(rotated_cyl)
    
    # Union all cylinders
    final_shape = geom.boolean_union(rotated_cylinders)
    
    mesh = geom.generate_mesh(dim=3)

import pygmsh
import meshio

# Create and optimize mesh
with pygmsh.geo.Geometry() as geom:
    circle = geom.add_circle([0, 0, 0], 1.0, mesh_size=0.1)
    mesh = geom.generate_mesh(dim=2)

# Optimize mesh quality
optimized_mesh = pygmsh.optimize(mesh, verbose=True)

# Write to multiple formats for different tools
pygmsh.write("mesh.msh")        # Gmsh format
pygmsh.write("mesh.vtk")        # VTK for ParaView  
pygmsh.write("mesh.med")        # MED for Code_Aster
pygmsh.write("mesh.stl")        # STL for 3D printing

# Extract to meshio for analysis
mesh_data = pygmsh.extract_to_meshio()

# Analyze mesh statistics
print(f"Number of nodes: {len(mesh_data.points)}")
print(f"Number of elements: {sum(len(cell.data) for cell in mesh_data.cells)}")
print(f"Physical groups: {list(mesh_data.cell_sets.keys())}")

# Compute mesh quality metrics
element_areas = []
for cell_block in mesh_data.cells:
    if cell_block.type == "triangle":
        # Compute triangle areas
        triangles = cell_block.data
        points = mesh_data.points
        for tri in triangles:
            p0, p1, p2 = points[tri]
            area = 0.5 * np.linalg.norm(np.cross(p1 - p0, p2 - p0))
            element_areas.append(area)

print(f"Average element area: {np.mean(element_areas):.6f}")
print(f"Element area std dev: {np.std(element_areas):.6f}")

import pygmsh
import numpy as np

def create_and_process_mesh(geometry_func, optimization_method="HighOrder"):
    """
    Complete mesh processing pipeline with geometry creation, 
    mesh generation, optimization, and quality analysis.
    """
    
    # Create geometry
    with pygmsh.occ.Geometry() as geom:
        # Apply user-defined geometry function
        geometry_entities = geometry_func(geom)
        
        # Generate initial mesh
        initial_mesh = geom.generate_mesh(dim=3)
    
    # Optimize mesh
    optimized_mesh = pygmsh.optimize(
        initial_mesh, 
        method=optimization_method, 
        verbose=True
    )
    
    # Extract mesh data for analysis
    mesh_data = pygmsh.extract_to_meshio()
    
    # Write mesh in multiple formats
    formats = ["mesh.msh", "mesh.vtk", "mesh.med"]
    for fmt in formats:
        pygmsh.write(fmt)
    
    return {
        'initial_mesh': initial_mesh,
        'optimized_mesh': optimized_mesh,
        'mesh_data': mesh_data,
        'num_nodes': len(mesh_data.points),
        'num_elements': sum(len(cell.data) for cell in mesh_data.cells)
    }

# Example usage
def complex_geometry(geom):
    """Create complex geometry for testing."""
    main_body = geom.add_ball([0, 0, 0], 2.0)
    
    holes = [
        geom.add_cylinder([x, y, -3], [0, 0, 1], 0.3)
        for x, y in [(-1, -1), (1, -1), (1, 1), (-1, 1)]
    ]
    
    return geom.boolean_difference(main_body, holes)

# Process mesh with pipeline
results = create_and_process_mesh(complex_geometry, "HighOrderElastic")
print(f"Generated mesh with {results['num_elements']} elements")

def optimize_cli(argv=None):
    """
    Command-line interface for mesh optimization.
    
    Provides a command-line tool for optimizing mesh files using
    various Gmsh optimization algorithms. Supports all meshio-compatible
    formats for input and output.
    
    Command line usage:
    pygmsh-optimize input.vtk output.xdmf [--method METHOD] [--verbose]
    
    Parameters:
    - argv: list[str], optional command line arguments (uses sys.argv if None)
    
    Supported formats:
    - Input: .msh, .vtk, .vtu, .stl, .ply, .obj, .off, .med
    - Output: .msh, .xdmf, .vtk, .vtu, .stl, .ply, .obj, .off, .med
    
    Returns:
    int: Exit code (0 for success, non-zero for errors)
    
    Examples:
    # Optimize VTK mesh and save as XDMF
    pygmsh-optimize mesh.vtk optimized.xdmf
    
    # Use specific optimization method with verbose output
    pygmsh-optimize input.msh output.msh --method "HighOrderOptimize" --verbose
    """