I am trying to do a boolean difference between two different shapes (a larger cylinder and a smaller cylinder that’s deformed).

The reasons is I want to use this geometry input for a simulation later on. Here’s my code:

import pyvista as pv
import numpy as np
import os
import json
from pymoab import core, types

height = 5
radius = 0.1
fuel_radius = 0.5  # Define the radius of the fuel region larger than the deformed cylinder
resolution_circumference = 40
num_parameters = 5
subdivs = 5
V_target = np.pi * radius ** 2 * height  # Arbitrary target volume for scaling
tol = height/100

def generate_deformed_cylinder(height, radius, resolution_circumference, num_parameters, V_target):
    cylinder = pv.Cylinder(center=(0.0, 0.0, 0.0), radius=radius, height=height,
                           resolution=resolution_circumference).triangulate().subdivide(subdivs)
    cylinder_polydata = cylinder.extract_surface().compute_normals(cell_normals=False)
    normals = cylinder_polydata['Normals']
    deformed_points = np.copy(cylinder_polydata.points)

    sigmas = np.random.uniform(0, height/2, size=num_parameters)
    sgns_scale = np.random.uniform(-radius/4, 0, size=num_parameters)
    centers = np.random.uniform(-height//2, height//2, size=num_parameters)

    for sgn, sigma, center in zip(sgns_scale, sigmas, centers):
        for i, point in enumerate(deformed_points):
            if point[0] > -height / 2 + tol and point[0] < height / 2 - tol:
                distance = np.linalg.norm(point - np.array([center, 0, 0]))
                gaussian_factor = np.exp(-0.5 * (distance / sigma) ** 2)
                displacement = sgn * gaussian_factor * normals[i]
                deformed_points[i] += displacement

    deformed_cylinder = pv.PolyData(deformed_points, faces=cylinder_polydata.faces)
    deformed_volume = deformed_cylinder.volume
    scaling_factor = (V_target / deformed_volume) ** (1 / 2)
    scaled_points = deformed_cylinder.points * [1, scaling_factor, scaling_factor]
    scaled_cylinder = pv.PolyData(scaled_points, deformed_cylinder.faces)

    return scaled_cylinder.triangulate()

def export_to_stl(mesh, filename):
    mesh.save(filename)
    print("Saved stl as {}".format(filename))

# Generate deformed control rod cylinder
control_rod = generate_deformed_cylinder(height, radius, resolution_circumference, num_parameters, V_target)
export_to_stl(control_rod, r"deformed_cylinder.stl")

# Create the larger fuel cylinder
fuel_cylinder = pv.Cylinder(center=(0.0, 0.0, 0.0), radius=fuel_radius, height=height, resolution=resolution_circumference).triangulate().subdivide(subdivs)
# Subtract the deformed control rod from the fuel cylinder
fuel_region = fuel_cylinder.extract_surface().boolean_difference(control_rod).clean()
export_to_stl(fuel_region, r"fuel_region.stl")

So I tried doing the boolean difference, and I exported this as stl files. I realized when the smaller cylinder is barely deformed, then I get a nice subtract of the geometries. However, when the cylinder becomes slightly deformed, then I get these results:

This s the result of the boolean_difference, and I get some weird artefacts on the top and bottom cap, but the final surface mesh isn’t even connected.

I am unsure what causes this issue, and I have tried looking for people with similiar issues. I found some threads about having difficulties with the boolean difference before, but it didn’t give the issues I see here.

Is there any obvious fault in my code?