Firstly, I am pleased with this new feature! Nevertheless, I did stubble upon an issue.
I am testing the example code provided on this page for passing mesh sizing functions as background meshes. The python script I am using is depicted below. When running this code with a sphere radius of 0.5 (default), everything works perfectly. When increasing the radius to 0.9 or bigger, I get a segmentation error:
_Initializing memorypools.
tetrahedron per block: 8188.
Size of a point: 168 bytes.
Size of a tetrahedron: 96 (96) bytes.
Size of a shellface: 200 (200) bytes.
Initializing robust predicates.
sizeof(double) = 8
machine epsilon = 2.22045e-16 [IEEE 754 64-bit macheps]
Delaunizing vertices...
Permuting vertices.
Sorting vertices.
Incrementally inserting vertices.
Delaunay seconds: 0.000397
Point sorting seconds: 1.4e-05
Creating surface mesh ...
160 (170) subfaces (segments).
Surface mesh seconds: 0.000192
Initializing memorypools.
tetrahedron per block: 8188.
Size of a point: 152 bytes.
Size of a tetrahedron: 96 (96) bytes.
Size of a shellface: 200 (200) bytes.
Initializing robust predicates.
sizeof(double) = 8
machine epsilon = 2.22045e-16 [IEEE 754 64-bit macheps]
Reconstructing mesh ...
Background mesh reconstruct seconds: 0.012948
Interpolating mesh size ...
Segmentation fault (core dumped)_
Code:
import pyvista as pv
import tetgen
import numpy as np
# Generate a background mesh with desired resolution
def generate_background_mesh(bounds, resolution=20, eps=1e-16):
x_min, x_max, y_min, y_max, z_min, z_max = bounds
grid_x, grid_y, grid_z = np.meshgrid(
np.linspace(x_min - eps, x_max + eps, resolution),
np.linspace(y_min - eps, y_max + eps, resolution),
np.linspace(z_min - eps, z_max + eps, resolution),
indexing="ij",)
return pv.StructuredGrid(grid_x, grid_y, grid_z).triangulate()
# Define sizing function based on proximity to a point of interest
def sizing_function(points, focus_point=np.array([0, 0, 0]), max_size=0.5, min_size=0.01):
distances = np.linalg.norm(points - focus_point, axis=1)
distances = np.clip(max_size - distances, min_size, max_size)
distances = np.ones(len(distances))*10
print(len(distances))
return distances
def mshVolumetricWriter(path, nodes, tets, labels):
# Open the .msh file for writing
with open(path, 'w') as file:
# Write the header
file.write("$MeshFormat\n2.2 0 8\n$EndMeshFormat\n")
# Write the nodes
file.write("$Nodes\n")
file.write(str(len(nodes)) + "\n")
for i, node in enumerate(nodes):
file.write("{} {:.6f} {:.6f} {:.6f}\n".format(i + 1, node[0], node[1], node[2]))
file.write("$EndNodes\n")
# Write the tetrahedra
file.write("$Elements\n")
file.write(str(len(tets)) + "\n")
for i, tet in enumerate(tets):
file.write("{} 4 2 0 {} {} {} {} {}\n".format(i + 1, labels[i], tet[0] + 1, tet[1] + 1, tet[2] + 1, tet[3] + 1))
file.write("$EndElements\n")
return 1
# Optionally write out the background mesh
def write_background_mesh(background_mesh, out_stem):
"""Write a background mesh to a file.
This writes the mesh in tetgen format (X.b.node, X.b.ele) and a X.b.mtr file
containing the target size for each node in the background mesh.
"""
mtr_content = [f"{background_mesh.n_points} 1"]
target_size = background_mesh.point_data["target_size"]
for i in range(background_mesh.n_points):
mtr_content.append(f"{target_size[i]:.8f}")
# Load or create your PLC
mesh = pv.Sphere(theta_resolution=10, phi_resolution=10, radius=10) # [CHANGE] radius of 0.5 works but larger radius does not
bg_mesh = generate_background_mesh(mesh.bounds)
bg_mesh.point_data['target_size'] = sizing_function(bg_mesh.points)
write_background_mesh(bg_mesh, 'bgmesh.b')
tet_kwargs = dict(verbose=True, order=1, mindihedral=20, minratio=1.5)
tet = tetgen.TetGen(mesh)
tet.tetrahedralize(bgmesh=bg_mesh, **tet_kwargs) # [ISSUE]
refined_mesh = tet.grid
nodes = np.array(refined_mesh.points)
tets = refined_mesh.cells.reshape(-1, 5)[:, 1:]
mshVolumetricWriter('/home/cjrodriguez/Sphere_refined.msh', nodes, tets, np.zeros(tets.shape[0], dtype=int))
Hi all,
Firstly, I am pleased with this new feature! Nevertheless, I did stubble upon an issue.
I am testing the example code provided on this page for passing mesh sizing functions as background meshes. The python script I am using is depicted below. When running this code with a sphere radius of 0.5 (default), everything works perfectly. When increasing the radius to 0.9 or bigger, I get a segmentation error:
Code: