When I generate a clump using a template, it always reports a volume of “1,” regardless of the clump’s size. For example, changing the radius at line 20 (radius = 2 # Radius of the sphere) does not affect the reported volume. As a result, the real mass always equals the density specified in the “clump generate” command.
I have attached an example Python script and a screenshot.
Regards,
get_ipython().magic('reset -sf') #command to initialize iPython environment
get_ipython().magic('clear') #iPython screen clear
import itasca as it
import os
it.command("python-reset-state false")
it.command("""
model new
model large-strain on
model domain extent -10 10 -10 10 -10 10
""")
#### Define Input variables ####
script_directory = r"D:\Clumps"
file_name = "clumps.txt"
path_to_file = os.path.join(script_directory, file_name)
radius = 2 # Radius of the sphere
num_latitude = 5 # Number of latitude bands
num_longitude = 5 # Number of longitude slices
randomness_factor = 0.1 # Factor by which to vary the radius randomly
bpDistance = 100 # 0 ~ 180 The greater, the smoother the pebble distribution
bpRatio = 0.3 # Ratio of the smallest to largest pebble kept in the clump template with 0 < fratio < 1.
bpRadFactor = 1.1 # to control pebble intrusion depth 1.05 is default > 1.0
bprefinenum = 10000 # Increase for dense pebble (Caution)
#######################################################
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D, art3d
def generate_sphere_mesh_with_randomness(radius, num_latitude, num_longitude, randomness_factor):
vertices = []
facets = []
# Top vertex
vertices.append((0, 0, radius))
# Generate vertices
for i in range(1, num_latitude):
theta = np.pi * i / num_latitude # Latitude angle
sin_theta = np.sin(theta)
cos_theta = np.cos(theta)
for j in range(num_longitude):
phi = 2 * np.pi * j / num_longitude # Longitude angle
rand_radius = radius * (1 + np.random.uniform(-randomness_factor, randomness_factor))
x = rand_radius * sin_theta * np.cos(phi)
y = rand_radius * sin_theta * np.sin(phi)
z = radius * cos_theta # z remains the same to maintain the overall shape
vertices.append((x, y, z))
# Bottom vertex
vertices.append((0, 0, -radius))
# Generate triangular facets for the top cap
for j in range(num_longitude):
p1 = 0
p2 = j + 1
p3 = (j + 1) % num_longitude + 1
facets.append((p1, p2, p3))
# Generate triangular facets for the middle
for i in range(1, num_latitude - 1):
for j in range(num_longitude):
p1 = (i - 1) * num_longitude + j + 1
p2 = p1 + num_longitude
p3 = p2 + 1 if (j + 1) < num_longitude else p2 + 1 - num_longitude
p4 = p1 + 1 if (j + 1) < num_longitude else p1 + 1 - num_longitude
facets.append((p1, p2, p3))
facets.append((p1, p3, p4))
# Generate triangular facets for the bottom cap
bottom_vertex_index = len(vertices) - 1
for j in range(num_longitude):
p1 = bottom_vertex_index
p2 = bottom_vertex_index - (j + 1)
p3 = bottom_vertex_index - ((j + 1) % num_longitude + 1)
facets.append((p1, p2, p3))
return vertices, facets
def compute_normal(v0, v1, v2):
edge1 = np.subtract(v1, v0)
edge2 = np.subtract(v2, v0)
normal = np.cross(edge1, edge2)
normal = normal / np.linalg.norm(normal) # Normalize the vector
return normal
def ensure_outward_normals(vertices, facets):
sphere_center = np.array([0, 0, 0])
corrected_facets = []
for facet in facets:
v0 = np.array(vertices[facet[0]])
v1 = np.array(vertices[facet[1]])
v2 = np.array(vertices[facet[2]])
normal = compute_normal(v0, v1, v2)
center = np.mean([v0, v1, v2], axis=0)
dot_product = np.dot(normal, center - sphere_center)
if dot_product < 0: # Normal is pointing inward
corrected_facets.append((facet[0], facet[2], facet[1]))
else:
corrected_facets.append(facet)
return corrected_facets
def verify_normals(vertices, facets):
normals = []
for facet in facets:
v0 = vertices[facet[0]]
v1 = vertices[facet[1]]
v2 = vertices[facet[2]]
normal = compute_normal(np.array(v0), np.array(v1), np.array(v2))
normals.append(normal)
center = np.mean([v0, v1, v2], axis=0)
dot_product = np.dot(normal, center)
print(f"Facet {facet} normal: {normal}, dot product with center: {dot_product}")
return normals
def plot_sphere_mesh_with_normals(vertices, facets, normals):
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Plot vertices
x = [v[0] for v in vertices]
y = [v[1] for v in vertices]
z = [v[2] for v in vertices]
ax.scatter(x, y, z, c='r', marker='o')
for i, (xi, yi, zi) in enumerate(vertices):
ax.text(xi, yi, zi, f'V{i}', color='blue')
# Plot normals
for i, normal in enumerate(normals):
facet = facets[i]
v0 = vertices[facet[0]]
v1 = vertices[facet[1]]
v2 = vertices[facet[2]]
center = np.mean([v0, v1, v2], axis=0)
ax.quiver(center[0], center[1], center[2], normal[0], normal[1], normal[2], length=1.0, color='blue')
# Plot facets
for facet in facets:
v0 = vertices[facet[0]]
v1 = vertices[facet[1]]
v2 = vertices[facet[2]]
poly3d = [[v0, v1, v2]]
ax.add_collection3d(art3d.Poly3DCollection(poly3d, facecolors='cyan', linewidths=1, edgecolors='r', alpha=.25))
plt.show()
def print_facet_vertices(vertices, facets):
for i, facet in enumerate(facets):
v0 = vertices[facet[0]]
v1 = vertices[facet[1]]
v2 = vertices[facet[2]]
print(f"Facet {i}:")
print(f" Vertex 0: {v0}")
print(f" Vertex 1: {v1}")
print(f" Vertex 2: {v2}")
# Generate the mesh
vertices, facets = generate_sphere_mesh_with_randomness(radius, num_latitude, num_longitude, randomness_factor)
# Ensure normals point outward
corrected_facets = ensure_outward_normals(vertices, facets)
# Verify the normals
normals = verify_normals(vertices, corrected_facets)
# Print the vertices of each facet
print_facet_vertices(vertices, corrected_facets)
# Plot the vertices, facets, and normals
#plot_sphere_mesh_with_normals(vertices, corrected_facets, normals)
for i, facet in enumerate(facets):
v0 = vertices[facet[0]]
v1 = vertices[facet[1]]
v2 = vertices[facet[2]]
it.command(f"""
geometry polygon create by-positions {v0} {v1} {v2} set 'whiba'
""")
it.command(f"""
clump template create name 'akas' ...
geometry 'whiba' ...
bubblepack ratio {bpRatio} distance {bpDistance} radfactor {bpRadFactor} refinenum {bprefinenum} ...
calculate-surface
clump generate density 2700 number 1 template 1 'akas'
""")
with open(path_to_file, 'w') as f:
f.write("# Molecule for Clumps\n")
f.write("# header section:\n")
f.write(f"{it.clump.pebble.count()} atoms\n")
clump_obj = it.clump.find(1)
f.write(f"{clump_obj.mass_real()} mass\n")
f.write(f"{clump_obj.pos_x()} {clump_obj.pos_y()} {clump_obj.pos_z()} com\n")
moi_real = clump_obj.moi_real()
f.write(f"{moi_real.xx()} {moi_real.yy()} {moi_real.zz()} {moi_real.xy()} {moi_real.xz()} {moi_real.yz()} inertia\n\n")
f.write(f"# body section:\n")
f.write(f"Coords\n\n")
pebble_list = it.clump.pebble.list() # Returns a list of pebble objects
for pebble_obj in pebble_list:
f.write(f"{pebble_obj.id()} {pebble_obj.pos_x()} {pebble_obj.pos_y()} {pebble_obj.pos_z()}\n")
f.write(f"\nDiameters\n\n")
pebble_list = it.clump.pebble.list() # Returns a list of pebble objects You must reset this!!!!
for pebble_obj in pebble_list:
diameter = 2.0 * pebble_obj.radius()
f.write(f"{pebble_obj.id()} {diameter}\n")
