diff --git a/apps/showcases/FreeSurface/raytracerVDB-1.py b/apps/showcases/FreeSurface/raytracerVDB-1.py
deleted file mode 100644
index a4cf0464e15c4debec39b20ae30963957cffa4e0..0000000000000000000000000000000000000000
--- a/apps/showcases/FreeSurface/raytracerVDB-1.py
+++ /dev/null
@@ -1,1019 +0,0 @@
-import pyopenvdb as vdb
-import os
-from math import floor
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib.colors as colors
-from matplotlib.animation import FuncAnimation
-from typing import Tuple, List
-import scipy.ndimage as ndimage
-from scipy.ndimage import gaussian_filter
-from scipy.spatial.transform import Rotation as Ro
-from scipy.spatial.transform import Slerp
-from scipy.interpolate import PchipInterpolator
-from sympy.printing.tree import print_node
-
-
-def runRayTracer(x_num, y_num, z_num, fill_levels, current_timestep):
-
- print("current_timestep", current_timestep)
- # fill_levels = np.clip(fill_levels, 0, 1)
- fill_levels = fill_levels.astype(np.float64)
-
- # Define the tolerance
- tolerance = 1e-8
-
- # Check if any value in fill_levels is outside the [0, 1] range
- out_of_range = np.logical_or(fill_levels < -tolerance, fill_levels > 1 + tolerance)
- if np.any(out_of_range):
- out_of_range_values = fill_levels[out_of_range]
- out_of_range_indices = np.where(out_of_range)
-
- for idx, value in zip(zip(*out_of_range_indices), out_of_range_values):
- print(f"Index: {idx}, Value: {value}")
-
- # raise ValueError(f"Fill level values must be between 0 and 1 inclusive (with tolerance of {tolerance}).")
-
- # Check for fill levels within the range [0, 1]
- in_range = np.logical_and(fill_levels > 0, fill_levels < 1)
- in_range_values = fill_levels[in_range]
- in_range_indices = np.where(in_range)
-
- # print("Fill level values within range (0 to 1):")
- # for idx, value in zip(zip(*in_range_indices), in_range_values):
- # if idx == (np.int64(309),):
- # print(f"Index: {idx}, Value: {value}")
-
- # Constants
- pi: float = np.pi # Mathematical constant π
- e: float = 1.60217662e-19 # Elementary charge in Coulombs
- beam_current: float = 3e-3 # Beam current in Amperes
- dwell_time: float = 0.32e-6 # Dwell time in seconds
- A: float = 21.5 # Approximation for Ti6Al4V
- # Ti6Al4V = Ti6Al4V.create_Ti6Al4V(T)
- # A: float = Ti6Al4V.atomic_number
- # print('Ti6Al4V.elements: ', Ti6Al4V.elements)
- # print('Ti6Al4V.composition: ', Ti6Al4V.composition)
- # print('A: ', A)
- dx: float = 3e-6 # Cell size in meters
- dy: float = dx # Cell size in meters
- dz: float = dx # Cell size in meters
- threshold: float = 0 # Beam profile threshold
- # The voltage at detector d can be computed using the transimpedance T (in units of V /A ) of a virtual amplifier in the context of the model
- transimpedance: float = 1e5 # Transimpedance in Volts per Ampere
- eV_threshold = 50 # Threshold energy in eV
- threshold_energy = eV_threshold * e # Threshold energy in joules
- # Define the work distance (height of the electron beam gun)
- # w_d = 12e-3 # Work distance in meters (272 mm)
-
- # Domain setup (Section 2.1.1)
- # Define the 3D domain for the simulation
- x_min, x_max = -0.5 * x_num * dx, 0.5 * x_num * dx # -10 mm to 10 mm
- # print('x_min, x_max:', x_min, x_max)
- y_min, y_max = -0.5 * y_num * dy, 0.5 * y_num * dy # -10 mm to 10 mm
- # print('y_min, y_max:', y_min, y_max)
- z_min, z_max = 0, z_num * dz
- # print('z_min, z_max:', z_min, z_max)
- # x_num = int((x_max - x_min) / dx) + 1
- # y_num = int((y_max - y_min) / dy) + 1
- # z_num = int((z_max - z_min) / dz) + 1
- # print(f"Domain dimensions: x={x_num}, y={y_num}, z={z_num}")
- # domain = np.zeros((z_num, y_num, x_num)) # (z, y , x)
- # print('domain shape: ', np.shape(domain))
- # domain[0:2, :, :] = 1 # Solid cells
- # domain[2, :, :] = 0.5 # Interface layer
- fill_levels_reshaped = np.reshape(fill_levels, (z_num, y_num, x_num))
- # fill_levels_reshaped = fill_levels
- # print('fill_levels_reshaped:', fill_levels_reshaped)
- out_of_range_fill_levels = np.logical_or(fill_levels_reshaped < -tolerance, fill_levels_reshaped > 1 + tolerance)
- if np.any(out_of_range_fill_levels):
- out_of_range_fill_level_values = fill_levels_reshaped[out_of_range_fill_levels]
- out_of_range_fill_level_indices = np.where(out_of_range_fill_levels)
-
- for idx, value in zip(zip(*out_of_range_fill_level_indices), out_of_range_fill_level_values):
- print(f"Index: {idx}, Value: {value}")
- raise ValueError(f"Fill level values must be between 0 and 1 inclusive (with tolerance of {tolerance}).")
-
- sigma = 0.0
- # Different sigma for each dimension
- # sigma = [0.0, 1.0, 0.5] # [z, y, x]
- fill_levels_smoothed = ndimage.gaussian_filter(fill_levels_reshaped, sigma=sigma, mode='nearest', radius=1)
-
- fill_levels_smoothed = np.clip(
- fill_levels_smoothed,
- a_min=np.float64(0.0),
- a_max=np.float64(1.0)
- )
-
- x = 242
- y = 189
- # Print fill levels at specific position
- print(f"\nFill levels at ({x}, {y}):")
- print("z_index | fill_level")
- print("--------------------")
- column = fill_levels_smoothed[:, y, x] # Get column at (157, 196)
- for z_idx, fill_level in enumerate(column):
- if 166 <= z_idx <= 176: # Print relevant range
- print(f" {z_idx} | {fill_level:.6f}")
-
- # Check if values are still in bounds
- if not (np.all(fill_levels_smoothed >= 0) and np.all(fill_levels_smoothed <= 1)):
- raise ValueError("Smoothed fill levels out of bounds!")
-
-
- def interpolate_normals(_normal_below: np.ndarray, _normal_above: np.ndarray, t: float) -> np.ndarray:
- """Interpolate between two normals using SLERP"""
- if not (0 <= t <= 1):
- raise ValueError(f"Interpolation parameter t must be in [0,1], got {t}")
-
- # Ensure normalized vectors
- _normal_below = normalize_vector(_normal_below.astype(np.float64))
- _normal_above = normalize_vector(_normal_above.astype(np.float64))
-
- # Compute dot product
- dot_product = np.clip(np.dot(_normal_below, _normal_above), -1.0, 1.0)
-
- # Thresholds for special cases
- PARALLEL_THRESHOLD = 0.9999
- ANTIPARALLEL_THRESHOLD = -0.9999
-
- # Case 1: Nearly parallel
- if dot_product > PARALLEL_THRESHOLD:
- return normalize_vector((1.0 - t) * _normal_below + t * _normal_above)
-
- # Case 2: Nearly antiparallel
- if dot_product < ANTIPARALLEL_THRESHOLD:
- perp = np.array([1.0, 0.0, 0.0])
- if abs(np.dot(_normal_below, perp)) > 0.9:
- perp = np.array([0.0, 1.0, 0.0])
- perp = normalize_vector(np.cross(_normal_below, perp))
- angle = np.pi * t
- return normalize_vector(_normal_below * np.cos(angle) + perp * np.sin(angle))
-
- # Case 3: Standard SLERP
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
-
- if abs(sin_theta) < 1e-10:
- return normalize_vector(_normal_below)
-
- scale0 = np.sin((1.0 - t) * theta) / sin_theta
- scale1 = np.sin(t * theta) / sin_theta
- return normalize_vector(scale0 * _normal_below + scale1 * _normal_above)
-
-
- def interpolate_fill_level_and_normal(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals: {surface_normals.ndim}")
-
- # Get fill levels at position
- _column = _fill_level[:, _y, _x]
-
- # Constants with physical meaning
- INTERFACE_THRESHOLD = 0.5 # Standard VOF interface threshold
- GRADIENT_THRESHOLD = 0.1 # Minimum gradient for interface detection
- TOLERANCE = 1e-6 # Numerical tolerance
-
- # Scan from top to bottom
- for z in range(len(_column)-2, -1, -1):
- current_fill = _column[z]
- next_fill = _column[z + 1]
-
- # Check for significant gradient
- gradient = abs(next_fill - current_fill)
- if gradient > GRADIENT_THRESHOLD:
- # Check for interface crossing (VOF method)
- '''crosses_interface = (
- (current_fill >= INTERFACE_THRESHOLD - TOLERANCE >= next_fill) or
- (current_fill <= INTERFACE_THRESHOLD + TOLERANCE <= next_fill)
- )'''
- crosses_interface = (
- (current_fill >= (INTERFACE_THRESHOLD - TOLERANCE) and
- next_fill <= (INTERFACE_THRESHOLD + TOLERANCE))
- )
-
- if crosses_interface:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}")
- print(f"fill_level_above: {next_fill}")
- # Linear interpolation with better numerical stability
- denominator = next_fill - current_fill
- if abs(denominator) < TOLERANCE:
- _z_interp = z + 0.5
- t_interface = 0.5
- else:
- # Calculate interpolation parameter
- t_interface = (INTERFACE_THRESHOLD - current_fill) / denominator
- if not (0 <= t_interface <= 1):
- error_msg = (
- f"\nInterpolation parameter t_interface={t_interface:.6f} is out of bounds [0,1]\n"
- f"Debug information:\n"
- f" - z-index: {z}\n"
- f" - current_fill: {current_fill:.6f}\n"
- f" - next_fill: {next_fill:.6f}\n"
- f" - denominator: {denominator:.6f}\n"
- f" - INTERFACE_THRESHOLD: {INTERFACE_THRESHOLD}\n"
- f" - Position (y,x): ({_y}, {_x})\n"
- f" - Fill level values around interface:\n"
- )
- # Nearby fill levels
- for i in range(max(0, z-2), min(len(_column), z+3)):
- error_msg += f" {i:3d} | {_column[i]:.6f}\n"
-
- raise ValueError(error_msg)
- # t_interface = np.clip(t_interface, 0.0, 1.0)
- _z_interp = z + t_interface
-
- print("_z_interp:", _z_interp)
- _z_height = (_z_interp + 0.5) * dz
- print("_z_height:", _z_height)
-
- # Get and validate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- _normal_interp = interpolate_normals(_normal_below, _normal_above, t_interface)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- return None, None, None, None, None, None, None
-
-
- def compute_surface_normals_sobel(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with default upward normal
- _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized[..., 2] = 1.0
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
-
- # Compute surface normals (Section 2.1.1)
- n_S = compute_surface_normals_sobel(fill_levels_smoothed)
- # print('np.shape(n_S): ', np.shape(n_S))
-
- # Beam profile (Section 2.1)
- # Define the Gaussian beam profile
- sigma: float = 500e-6 / (4 * dx) # Beam profile sigma
- x = np.arange(-np.floor(min(x_num//4, 2*sigma)), 1+np.floor(min((x_num//4), 2*sigma)))
- # print('x:', x)
- y = np.arange(-floor(min(y_num//4, 2*sigma)), 1+floor(min((y_num//4), 2*sigma)))
- # print('y:', y)
- X, Y = np.meshgrid(x, y)
- beam_profile = np.exp(-(X ** 2 + Y ** 2) / (2 * sigma ** 2))
- beam_profile /= np.sum(beam_profile) # Normalize beam profile
- print('beam_profile shape: ', np.shape(beam_profile))
-
- beam_size = np.shape(beam_profile)[0] * np.shape(beam_profile)[1]
- # print('Beam size:', beam_size)
-
- # Detector setup (Section 2.1.1)
- # Define positions of the detectors
- detector_height = 4 * z_max # Detector height in meters (272 mm)
- detector_spacing_x = 1.5 * dx * x_num # Detector spacing in meters (127 mm)
- detector_spacing_y = detector_spacing_x # 1.5 * dy * y_num # Detector spacing in meters (127 mm)
- detector_diameter = 4 * dx # Detector diameter in meters (50 mm)
-
- # Define detector positions relative to the origin
- detector_positions = np.array([
- [0.5 * detector_spacing_x, 0, detector_height],
- [-0.5 * detector_spacing_x, 0, detector_height],
- [0, 0.5 * detector_spacing_y, detector_height],
- [0, -0.5 * detector_spacing_y, detector_height]
- ], dtype=np.float64)
- # print('detector_positions:', detector_positions)
-
-
- # Define detector vertices for solid angle calculation
- def get_detector_vertices(detector_position: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
- # Calculate the side length of the detector
- side_length : float = np.sqrt(3) * (detector_diameter / 2)
- height : float = np.sqrt(3) * (side_length / 2)
-
- # Determine the orientation based on the detector position
- if np.allclose(detector_position[:2], [0, 0]):
- raise ValueError("Detector position cannot be at the origin (0, 0, 0)")
- elif detector_position[0] > 0: # Right detector
- _v1 = detector_position + np.array([2 / 3 * height, 0, 0], dtype=np.float64)
- _v2 = detector_position + np.array([-height / 3, side_length / 2, 0], dtype=np.float64)
- _v3 = detector_position + np.array([-height / 3, -side_length / 2, 0], dtype=np.float64)
- elif detector_position[0] < 0: # Left detector
- _v1 = detector_position + np.array([-2 / 3 * height, 0, 0], dtype=np.float64)
- _v2 = detector_position + np.array([height / 3, -side_length / 2, 0], dtype=np.float64)
- _v3 = detector_position + np.array([height / 3, side_length / 2, 0], dtype=np.float64)
- elif detector_position[1] > 0: # Back detector
- _v1 = detector_position + np.array([0, 2 / 3 * height, 0], dtype=np.float64)
- _v2 = detector_position + np.array([-side_length / 2, -height / 3, 0], dtype=np.float64)
- _v3 = detector_position + np.array([side_length / 2, -height / 3, 0], dtype=np.float64)
- elif detector_position[1] < 0: # Front detector
- _v1 = detector_position + np.array([0, -2 / 3 * height, 0], dtype=np.float64)
- _v2 = detector_position + np.array([side_length / 2, height / 3, 0], dtype=np.float64)
- _v3 = detector_position + np.array([-side_length / 2, height / 3, 0], dtype=np.float64)
- else:
- raise ValueError("Invalid detector position. Must be in one of the four expected positions.")
-
- return _v1, _v2, _v3
-
-
- # Compute solid angle using Equation A1 (Appendix A)
- def compute_solid_angle(vector1: np.ndarray, vector2: np.ndarray, vector3: np.ndarray) -> float:
- # Input validation
- if not isinstance(vector1, np.ndarray) or not isinstance(vector2, np.ndarray) or not isinstance(vector3, np.ndarray):
- raise ValueError("Inputs must be numpy arrays")
- if vector1.shape != (3,) or vector2.shape != (3,) or vector3.shape != (3,):
- raise ValueError("Inputs must be 3D vectors")
- if not np.isfinite(vector1).all() or not np.isfinite(vector2).all() or not np.isfinite(vector3).all():
- raise ValueError("Inputs must be finite")
-
- # Convert to high precision once
- vectors = [vec.astype(np.float64) for vec in (vector1, vector2, vector3)]
- cross_product = np.cross(vectors[1], vectors[2])
- # Check if the cross product is nearly zero
- if np.linalg.norm(cross_product) < 1e-16:
- raise ValueError("Vectors are nearly coplanar")
-
- # Compute solid angle
- try:
- # Compute the numerator of the solid angle formula
- numerator: float = np.dot(vectors[0], cross_product)
- # Check if the numerator is nearly zero
- if np.abs(numerator) < 1e-16:
- raise ValueError("Vectors are nearly parallel")
-
- norms = [np.linalg.norm(v) for v in vectors]
-
- # Compute the denominator of the solid angle formula
- denominator = (
- norms[0] * norms[1] * norms[2] +
- np.dot(vectors[0], vectors[1]) * norms[2] +
- np.dot(vectors[0], vectors[2]) * norms[1] +
- np.dot(vectors[1], vectors[2]) * norms[0]
- )
-
- # Compute the solid angle in radians
- solid_angle = 2 * np.arctan2(numerator, denominator)
- # print("Solid angle:", solid_angle)
- return solid_angle
- except ZeroDivisionError:
- raise ValueError("Error computing solid angle: division by zero")
- except np.linalg.LinAlgError:
- raise ValueError("Error computing solid angle: linear algebra error")
-
-
- # BSE coefficient functions (Equations 11 and 10, Section 2.1.2)
- def eta_0(_A: float) -> float:
- """Equation 11 from the paper"""
- _eta_0: float = -0.0254 + 0.016 * _A - 1.86e-4 * _A ** 2 + 8.3e-7 * _A ** 3
- print('eta_0: ', _eta_0)
- return _eta_0
-
-
- # BSE Coefficient for single atomic targets, found to hold in the range of 10–100 keV
- def eta(theta: float, _A: float) -> float:
- """Equation 10 from the paper"""
- B: float = 0.89
- _eta: float = B * (eta_0(_A) / B) ** np.cos(theta)
- print('eta: ', _eta)
- return _eta
-
-
- '''
- # Compute the weighted average BSE coefficient for the alloy
- def compute_weighted_eta(theta: float, alloy: Alloy) -> float:
- """Equation 12 from the paper"""
- weighted_eta: float = 0
- for element, weight_fraction in zip(alloy.elements, alloy.composition):
- # print('element: ', element)
- # print('weight_fraction: ', weight_fraction)
- _A: float = element.atomic_number
- # print('A: ', _A)
- eta_i: float = eta(theta, _A)
- # print('eta_i: ', eta_i)
- weighted_eta += weight_fraction * eta_i
- # print('weighted_eta: ', weighted_eta)
- return weighted_eta
- '''
-
-
- # Correction function C(theta) (Equation 14, Section 2.1.3)
- def C(theta: float) -> float:
- print("theta1:", theta)
- """Equation 14 from the paper"""
- theta = np.rad2deg(theta) # Ensure angle is in degrees for C function
- _C: float = 1.01 - 4.06 * theta + 1.36e-5 * theta ** 2 - 1.71e-6 * theta ** 3
- print('C: ', _C)
- # Check if C is negative and raise an exception
- if _C < 0:
- raise ValueError(f"C function returned a negative value for theta: {theta}. "
- f"Please check the input and the C function implementation.")
- return _C
-
-
- # Scattering function g_BSE (Equation 13, Section 2.1.3)
- def g_BSE(_theta: float, _alpha: float, _beta: float) -> float:
- print("_theta, _alpha, _beta:", _theta, _alpha, _beta)
- """Equation 13 from the paper"""
- # k: float = np.rad2deg(_theta) / 13 # Correct implementation as per the paper
- # print("k: ", k)
- diffusive_part: float = (eta_0(A) / pi) * np.cos(_alpha)
- print(np.cos(_alpha))
- print("diffusive_part:", diffusive_part)
- # reflective_part: float = ((eta(_theta, A) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
- # if np.abs(reflective_part) < 1e-10:
- # reflective_part = 0
- # print("reflective_part:", reflective_part)
-
- # Compute a weighted average BSE coefficient for the alloy
- # reflective_part: float = ((compute_weighted_eta(_theta, Ti6Al4V) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
-
- # print('g_BSE: ', diffusive_part + reflective_part)
- return diffusive_part #+ reflective_part
-
-
- # Calculate the beam direction vector p_E
- def calculate_p_E(_beam_x, _beam_y):
- print(f"calculate_p_E({_beam_x}, {_beam_y}):")
- # Calculate the vector from the beam gun to the beam location
- beam_vector = np.array([_beam_x * dx, _beam_y * dy, -detector_height], dtype=np.float64)
- return normalize_vector(beam_vector)
-
- def calculate_p_E1(_beam_x, _beam_y, _height):
- """Calculate normalized beam direction vector from gun to surface point
-
- Args:
- _beam_x: Beam x-position relative to center (in cells)
- _beam_y: Beam y-position relative to center (in cells)
- _height: Height difference between gun and surface point (m)
- """
- print(f"calculate_p_E({_beam_x}, {_beam_y}):")
- print("_height:", _height)
-
- # Calculate vector components
- x_component = _beam_x * dx # Convert cells to meters
- y_component = _beam_y * dy # Convert cells to meters
- z_component = _height # Already in meters
-
- # Create beam vector (from gun to surface point)
- beam_vector = np.array([x_component, y_component, z_component], dtype=np.float64)
-
- # Verify downward direction
- if normalize_vector(beam_vector)[2] > 0:
- raise ValueError("Beam vector pointing upward")
-
- # Normalize and return
- return normalize_vector(beam_vector)
-
- def calculate_p_E2(_beam_x, _beam_y, _height):
- """Calculate normalized beam direction vector from gun to surface point
-
- Args:
- _beam_x: x-position relative to center (in cells)
- _beam_y: y-position relative to center (in cells)
- _height: height difference between gun and surface (in meters)
-
- Returns:
- np.ndarray: Normalized beam direction vector
- """
- print(f"calculate_p_E({_beam_x}, {_beam_y}):")
- print("_height:", _height)
-
- # Calculate beam vector components (from gun to surface point)
- beam_vector = np.array([
- _beam_x * dx, # x component
- _beam_y * dy, # y component
- _height # z component (already in meters)
- ], dtype=np.float64)
-
- # Normalize the vector
- norm = np.linalg.norm(beam_vector)
- if norm < 1e-10:
- raise ValueError("Beam vector magnitude too small")
-
- normalized_vector = beam_vector / norm
-
- # Verify the vector is pointing downward (negative z)
- if normalized_vector[2] > 0:
- raise ValueError("Beam vector must point downward (negative z)")
-
- return normalized_vector
-
-
- # Calculate phi as the angle between p_E and the negative z-axis
- def calculate_phi(_p_E) -> float:
- z_axis = np.array([0.0, 0.0, -1.0], dtype=np.float64)
-
- # Check if _p_E is approximately equal to the negative z-axis
- '''if np.allclose(_p_E, z_axis, rtol=1e-8, atol=1e-8):
- _phi_rad = 0.0
- else:
- # Calculate the angle between _p_E and the z-axis using arccos
- _phi_rad = calculate_angle(_p_E, z_axis)'''
- _phi_rad = calculate_angle(_p_E, z_axis)
- print(f"calculate_phi({_p_E}, {_phi_rad}):")
- return _phi_rad
-
-
- def set_beam_movement(_x_center, _y_center, _x_offset, _y_offset, _movement_type):
- if _movement_type == 'X':
- if _x_offset >= 0:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset - 1
- else:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset + 1
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center + _y_offset + 1
- elif _movement_type == 'Y':
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center + _x_offset + 1
- if _y_offset >= 0:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset - 1
- else:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset + 1
- elif _movement_type == 'XY':
- if _x_offset >= 0:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset - 1
- else:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset + 1
- if _y_offset >= 0:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset - 1
- else:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset + 1
- else:
- raise ValueError("movement_type must be 'X' or 'Y' or 'XY'")
-
- return _beam_x_start, _beam_x_end, _beam_y_start, _beam_y_end
-
-
- # Initialize arrays to store electron counts
- electrons_per_detector_per_step = np.zeros((x_num, y_num, len(detector_positions)))
- total_electrons_per_step = np.zeros((x_num, y_num))
-
- # Additional arrays to capture detailed information
- electrons_per_detector_per_beam_profile = np.zeros(
- (beam_profile.shape[0], beam_profile.shape[1], len(detector_positions)))
- total_electrons_per_beam_profile = np.zeros((beam_profile.shape[0], beam_profile.shape[1]))
-
- # Calculate the center indices
- x_center = x_num // 2
- y_center = y_num // 2
- # print('x_center, y_center:', x_center, y_center)
-
- x_offset: int = 0
- y_offset: int = 0
- movement_type: str = 'Y'
-
- # Beam X-Movement
- # beam_x_start, beam_x_end = x_center - x_offset, x_center + x_offset + 1
- # beam_y_start, beam_y_end = y_center + y_offset, y_center + y_offset + 1
-
- # Beam Y-Movement
- # beam_x_start, beam_x_end = x_center + x_offset, x_center + x_offset + 1
- # beam_y_start, beam_y_end = y_center - y_offset, y_center + y_offset + 1
-
- # For X or Y Movement
- beam_x_start, beam_x_end, beam_y_start, beam_y_end = set_beam_movement(
- x_center, y_center, x_offset, y_offset, movement_type)
-
- # Determine the step size based on the direction of movement
- step_x = 1 if beam_x_start <= beam_x_end else -1
- step_y = 1 if beam_y_start <= beam_y_end else -1
-
- # print(beam_x_start, beam_x_end, beam_y_start, beam_y_end, step_x, step_y)
-
-
- z_heights = np.zeros(beam_size)
- fill_level_b = np.zeros(beam_size)
- fill_level_a = np.zeros(beam_size)
- normals = np.zeros((beam_size, 3))
- normals_b = np.zeros((beam_size, 3))
- normals_a = np.zeros((beam_size, 3))
- cell_counter = 0
-
-
- def normalize_vector(vector: np.ndarray) -> np.ndarray:
- """Normalize a 3D vector."""
- vector = vector.astype(np.float64) # Convert to float64
- norm = np.linalg.norm(vector)
- if norm < 1e-18: # Avoid division by zero
- return vector # Return original if nearly zero
- return vector / norm # Normalize
-
-
- def calculate_angle(_v1: np.ndarray, _v2: np.ndarray) -> float:
- """Calculate the angle between two normalized vectors."""
- _v1 = normalize_vector(_v1)
- _v2 = normalize_vector(_v2)
-
- dot_product = np.clip(np.dot(_v1, _v2), -1.0, 1.0)
-
- return np.arccos(dot_product) # Return angle in radians
-
-
- # Main simulation loop to iterate over the simulation domain
- for beam_x in range(beam_x_start, beam_x_end, step_x):
- for beam_y in range(beam_y_start, beam_y_end, step_y):
- beam_loc = (beam_x, beam_y)
- print(f"Beam center is at domain location: {beam_loc}")
-
- # Calculate the beam direction vector p_E
- # >> p_E = calculate_p_E(beam_x - x_center, beam_y - y_center)
- # >> print('p_E: ', p_E)
-
- # Calculate phi as the angle between p_E and the negative z-axis
- # phi_rad = np.arctan2(np.sqrt(dx**2 + dy**2), w_d)
- # >> phi_rad = calculate_phi(p_E)
- # print(f"phi_rad: {phi_rad}")
-
- # Reset counters for this beam location
- total_electrons: int = 0
- electrons_per_detector = np.zeros(len(detector_positions))
-
- # Iterate over the beam profile
- for profile_x in [125]:
- # for profile_x in range(beam_profile.shape[0]):
- # for profile_x in [beam_profile.shape[0]//2 + 0]:
- for profile_y in [72]:
- # for profile_y in range(beam_profile.shape[1]):
- # for profile_y in [beam_profile.shape[1]//2 + 0]:
- print(f"Beam profile shape ({profile_x}, {profile_y})")
- domain_x = beam_loc[0] - beam_profile.shape[0] // 2 + profile_x
- domain_y = beam_loc[1] - beam_profile.shape[1] // 2 + profile_y
- print(f"Domain ({domain_x}, {domain_y})")
-
- # p_E = calculate_p_E(domain_x - x_center, domain_y - y_center)
- # print('p_E: ', p_E)
- # phi_rad = calculate_phi(p_E)
-
- if (0 <= domain_x < x_num) and (0 <= domain_y < y_num):
- if beam_profile[profile_x, profile_y] > threshold:
- print('finding interface_z_indices')
- # print('interface_z_indices:', interface_z)
- z_interp, z_height, fill_level_below, fill_level_above, normal_interp, normals_below, normals_above \
- = interpolate_fill_level_and_normal(fill_levels_smoothed, n_S, domain_y, domain_x)
- if z_interp is not None:
- z_heights[cell_counter] = z_height
- fill_level_b[cell_counter] = fill_level_below
- fill_level_a[cell_counter] = fill_level_above
- normals[cell_counter] = normal_interp
- normals_b[cell_counter] = normals_below
- normals_a[cell_counter] = normals_above
- cell_counter += 1
- n_S_local = normal_interp
- print('n_S_local:', n_S_local)
- if not np.allclose(np.linalg.norm(n_S_local), 1.0, atol=1e-5):
- raise ValueError(f"n_S_local is not normalized, {np.linalg.norm(n_S_local)}")
-
- print("detector_height:", detector_height)
- print("z_height:", z_height)
- p_E = calculate_p_E1(domain_x - x_center, domain_y - y_center, z_height-detector_height)
- print('p_E: ', p_E)
-
- theta_local = calculate_angle(-p_E, n_S_local)
- print('theta_local:', theta_local)
- # if theta_local < 0.0 or theta_local > 0.011:
- # raise ValueError(f"theta_local value should be 0.0, found {theta_local}")
- # theta_local: 0.016
-
- scattering_point = np.array([
- (domain_x - beam_x) * dx,
- (domain_y - beam_y) * dy,
- z_height
- ], dtype=np.float64)
- print(f'scattering_point ({domain_x}, {domain_y}, {z_interp}):', scattering_point)
-
- for det_idx, det_pos in enumerate(detector_positions):
- # print(f"Detector {det_idx + 1} position: {det_pos}")
- det_idx: int
- # det_pos: np.ndarray[float]
- print(f'det_pos[{det_idx + 1}]: {det_pos}')
-
- # Calculate the vector from the scattering point to the detector
- det_vec = det_pos - scattering_point
- det_vec = normalize_vector(det_vec)
- print(f'det_vec[{det_idx + 1}]: {det_vec}')
- if not np.allclose(np.linalg.norm(det_vec), 1.0):
- raise ValueError("det_vec is not normalized.")
-
- # alpha = np.arccos(np.clip(np.dot(det_vec, n_S_local), -1.0, 1.0))
- alpha = calculate_angle(det_vec, n_S_local)
- print(f'alpha: {alpha}')
-
- if alpha > np.pi / 2:
- continue
-
- p_E_reflected = p_E - 2.0 * np.dot(p_E, n_S_local) * n_S_local
- # p_E_reflected = np.where(np.abs(p_E_reflected) < 1e-8, 0.0, p_E_reflected)
- p_E_reflected = normalize_vector(p_E_reflected)
- print(f'p_E_reflected: {p_E_reflected}')
-
- # beta = np.arccos(np.clip(np.dot(p_E_reflected, det_vec), -1.0, 1.0))
- beta = calculate_angle(p_E_reflected, det_vec)
- print(f'beta: {beta}')
-
- # Apply energy threshold
- # electron_energy = N_0 * e * g_bse_value
- # print('electron_energy: ', electron_energy)
-
- # Calculate the initial energy of the primary electrons (E_0)
- V_acc = 60e3 # Accelerating voltage in Volts
- E_0 = V_acc * e # Initial energy of primary electrons in joules
-
- # Calculate the energy of the backscatter electrons
- E_BSE = E_0 * eta(theta_local, A)
- # Compute a weighted average BSE coefficient for the alloy
- # E_BSE = E_0 * compute_weighted_eta(theta_local, Ti6Al4V)
- # print('E_BSE: ', E_BSE)
-
- # Bifurcation logic for PE and SE
- # if E_BSE > threshold_energy:
- # Backscatter electrons (BSE)
- # Compute scattering function g_bse_value
- g_bse_value = g_BSE(theta_local, alpha, beta)
- print(f'g_bse_value: {g_bse_value}')
-
- # Ensure g_bse_value is non-negative
- if g_bse_value <= 0:
- raise ValueError(f"g_bse_value is {g_bse_value} (zero negative), check calculations.")
-
- # Compute solid angle (Equation A1)
- v1, v2, v3 = get_detector_vertices(det_pos)
- normal = np.cross(v2 - v1, v3 - v1)
- normal = normalize_vector(normal)
-
- dOmega = compute_solid_angle(v1 - scattering_point,
- v2 - scattering_point,
- v3 - scattering_point)
- print(f"Detector {det_idx + 1} solid angle: {dOmega}")
-
- # Compute number of electrons collected by the detector (Equation A2)
- # The beam ray stays on each grid point ij for the dwell time t_d and distributes N_0 PE into it
- N_0 = beam_current * dwell_time / e # Equation 3, Section 2.1.1
- print("N_0:", N_0)
-
- # There are N1 electrons in the chamber after the 1st scattering event
- N_1 = N_0 * eta(theta_local, A) # Equation 4, Section 2.1.1
- # Compute a weighted average BSE coefficient for the alloy
- # N_1 = N_0 * compute_weighted_eta(theta_local, Ti6Al4V) # Equation 4, Section 2.1.1
- print("N_1:", N_1)
-
- # The number of electrons N_1_d collected by a general detector d
- N_1_d = N_1 * np.sum(g_bse_value * dOmega) # Equation 5, Section 2.1.1
-
- # Accumulate electrons per detector for this beam profile point
- electrons_per_detector_per_beam_profile[profile_x, profile_y, det_idx] += N_1_d
-
- # Accumulate electrons per detector for this domain point
- electrons_per_detector_per_step[beam_x, beam_y, det_idx] += N_1_d
-
- # Add to total electrons for this beam profile point
- total_electrons_per_beam_profile[profile_x, profile_y] += N_1_d
-
- # Add to total electrons for this domain point
- total_electrons_per_step[beam_x, beam_y] += N_1_d
-
- # Accumulate electrons per detector
- electrons_per_detector[det_idx] += N_1_d
-
- # Add to total electrons
- total_electrons += N_1_d
-
- # print(f"After addition: electrons_per_detector[{det_idx + 1}]: {electrons_per_detector[det_idx]}")
- else:
- print(f"No interface found at ({domain_x}, {domain_y})")
-
- # Print detailed information for this beam profile point
- # print(f"Electrons per detector at this beam profile point ({profile_x}, {profile_y}): {electrons_per_detector_per_beam_profile[profile_x, profile_y]}")
- # print(f"Total electrons at this beam profile point ({profile_x}, {profile_y}): {total_electrons_per_beam_profile[profile_x, profile_y]}")
-
- # Print detailed information for this domain point after iterating over the beam profile
- # print(f"Electrons per detector at this domain point ({beam_x}, {beam_y}): {electrons_per_detector_per_step[beam_x, beam_y]}")
- # print(f"Total electrons at this domain point ({beam_x}, {beam_y}): {total_electrons_per_step[beam_x, beam_y]}")
-
- # print("Shape of electrons_per_detector_per_step:", electrons_per_detector_per_step.shape)
- # print("Number of lines:", len(lines))
-
- # Compute final results
- total_electrons_all_steps = np.sum(total_electrons_per_step)
- total_electrons_per_detector = np.sum(electrons_per_detector_per_step, axis=(0, 1))
-
- # Print final results
- # print(f"Total electrons hitting all detectors: {total_electrons_all_steps}")
- # for i, electrons in enumerate(total_electrons_per_detector):
- # print(f"Total electrons hitting detector {i + 1}: {electrons}")
-
- # for beam_x in range(beam_x_start, beam_x_end, step_x):
- # for beam_y in range(beam_y_start, beam_y_end, step_y):
- # print(f"Electrons per detector at this domain point ({beam_x}, {beam_y}): {electrons_per_detector_per_step[beam_x, beam_y]}")
- # print('end simulation')
- # print(electrons_per_detector_per_step)
-
- # print("total_electrons:", total_electrons, "at timestep", current_timestep)
- # print("electrons_per_det:", electrons_per_det, "at timestep", current_timestep)
- # print("electrons_per_beam_loc:", electrons_per_beam_loc, "at timestep", current_timestep)
- # print("electrons_per_beam_loc_per_det:", electrons_per_beam_loc_per_det, "at timestep", current_timestep)
- return electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, normals, normals_b, normals_a
-
-def read_fill_levels(filepath):
- """Read only the FillingFrac grid from VDB file"""
- try:
- # Try to read just the FillingFrac grid
- grid = vdb.read(filepath, "FillingFrac")
-
- # Get dimensions
- dims = grid.evalActiveVoxelDim()
- print(f"Grid dimensions: {dims}")
-
- # Create array
- fill_array = np.zeros(dims, dtype=np.float64)
-
- # Copy data
- grid.copyToArray(fill_array)
-
- # Transpose to ray tracer's expected order (z, y, x)
- fill_array = np.transpose(fill_array, (2, 1, 0))
-
- # Create standardized array with 256 z cells
- standardized_array = np.zeros((256, dims[1], dims[0]), dtype=np.float64)
-
- # Copy data from original array to standardized array
- standardized_array[:fill_array.shape[0], :, :] = fill_array
-
- # Apply Gaussian smoothing
- sigma = 1.0
- standardized_array_smoothed = ndimage.gaussian_filter(standardized_array, sigma=sigma, mode='nearest', radius=1)
-
- # Print values
- x = 242
- y = 189
- print(f"\nFill levels at ({x}, {y}):")
- print("z_index | fill_level")
- print("-" * 20)
-
- # Get fill levels
- fill_levels = standardized_array_smoothed[:, y, x]
- valid_levels = fill_levels[(fill_levels > 0) & (fill_levels < 1)]
-
- if len(valid_levels) > 0:
- for z, val in enumerate(fill_levels):
- if 0 <= val <= 1:
- print(f"{z:7d} | {val:.6f}")
- else:
- print("No fill level found between 0 and 1")
-
- return standardized_array_smoothed
-
- except Exception as e:
- print(f"Error reading {filepath}")
- print(f"Error details: {str(e)}")
- try:
- metadata = vdb.readMetadata(filepath)
- print("\nFile metadata:")
- for key, value in metadata.items():
- print(f"{key}: {value}")
- except Exception as meta_e:
- print(f"Error reading metadata: {str(meta_e)}")
- return None
-
-def process_single_timestep(filepath, timestep):
- """Process a single VDB file and run ray tracer"""
- try:
- # Read the FillingFrac grid
- grid = vdb.read(filepath, "FillingFrac")
-
- # Get dimensions
- dims = grid.evalActiveVoxelDim()
- x_num, y_num = dims[0], dims[1]
- z_num = 256 # dims[2] # Fixed z dimension
- print(f"Domain dimensions: {x_num} x {y_num} x {z_num}")
-
- # Read fill levels
- fill_levels = read_fill_levels(filepath)
- if fill_levels is None:
- raise ValueError("Could not read fill levels from VDB file")
-
- if fill_levels is not None:
- print("\nSuccessfully read fill levels")
- print(f"Array shape: {np.shape(fill_levels)}")
-
- try:
- # Convert to numpy array if not already
- fill_array = np.asarray(fill_levels)
-
- # Get non-zero elements
- non_zero = fill_array[fill_array > 0]
-
- # Get interface cells (0 < fill_level < 1)
- interface_cells = fill_array[(fill_array > 0) & (fill_array < 1)]
-
- # Print statistics
- if len(non_zero) > 0:
- print("\nNon-zero cells statistics:")
- print(f"Number of non-zero cells: {len(non_zero)}")
- print(f"Min non-zero value: {np.min(non_zero):.6f}")
- print(f"Max value: {np.max(fill_array):.6f}")
-
- if len(interface_cells) > 0:
- print("\nInterface cells statistics:")
- interface_percentage = (len(interface_cells) * 100.0) / len(non_zero)
- print(f"Number of interface cells: {len(interface_cells)} ({interface_percentage:.3f}%)")
- print(f"Min interface value: {np.min(interface_cells):.6f}")
- print(f"Max interface value: {np.max(interface_cells):.6f}")
- print(f"Mean interface value: {np.mean(interface_cells):.6f}")
- else:
- print("\nNo interface cells found (no fill levels between 0 and 1)")
-
- except Exception as e:
- print(f"Error processing array: {str(e)}")
-
- # Run ray tracer
- print(f"\nProcessing timestep {timestep}")
- electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, \
- normals, normals_b, normals_a = runRayTracer(
- x_num,
- y_num,
- z_num,
- # fill_levels.flatten(),
- fill_levels,
- timestep
- )
-
- return electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, \
- normals, normals_b, normals_a
-
- except Exception as e:
- print(f"Error processing timestep {timestep}: {str(e)}")
- raise
-
-def main():
- """Main function to process single VDB file and run ray tracer"""
- # Specific VDB file path
- filepath = "/local/ca00xebo/softwares/KiSSAM/Markl-plate/geometry-000006200.vdb"
- timestep = 0 # or extract from filename if needed
-
- try:
- print(f"\nProcessing file: {os.path.basename(filepath)}")
-
- # Process the file
- result = process_single_timestep(filepath, timestep)
-
- # Extract and process results
- electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, \
- normals, normals_b, normals_a = result
-
- # Calculate total electrons for each detector
- total_electrons_for_detectors = np.sum(electrons_per_detector_per_step, axis=(0, 1))
- print("Total electrons per detector:", total_electrons_for_detectors)
-
- return result
-
- except Exception as e:
- print(f"Error processing file: {str(e)}")
- raise
-
-
-if __name__ == "__main__":
- results = main()
\ No newline at end of file
diff --git a/apps/showcases/FreeSurface/raytracerVDB-120.py b/apps/showcases/FreeSurface/raytracerVDB-120.py
deleted file mode 100644
index 29ec1af8361d97f52b7d5e35a621cccc4ad2ff42..0000000000000000000000000000000000000000
--- a/apps/showcases/FreeSurface/raytracerVDB-120.py
+++ /dev/null
@@ -1,1737 +0,0 @@
-import pyopenvdb as vdb
-import os
-from math import floor
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib.colors as colors
-from matplotlib.animation import FuncAnimation
-from typing import Tuple, List
-import scipy.ndimage as ndimage
-from scipy.ndimage import gaussian_filter
-from scipy.spatial.transform import Rotation as Ro
-from scipy.spatial.transform import Slerp
-from scipy.interpolate import PchipInterpolator
-from sympy.printing.tree import print_node
-
-
-def runRayTracer(x_num, y_num, z_num, fill_levels, current_timestep):
-
- print("current_timestep", current_timestep)
- # fill_levels = np.clip(fill_levels, 0, 1)
- fill_levels = fill_levels.astype(np.float64)
-
- # Define the tolerance
- tolerance = 1e-8
-
- # Check if any value in fill_levels is outside the [0, 1] range
- out_of_range = np.logical_or(fill_levels < -tolerance, fill_levels > 1 + tolerance)
- if np.any(out_of_range):
- out_of_range_values = fill_levels[out_of_range]
- out_of_range_indices = np.where(out_of_range)
-
- for idx, value in zip(zip(*out_of_range_indices), out_of_range_values):
- print(f"Index: {idx}, Value: {value}")
-
- # raise ValueError(f"Fill level values must be between 0 and 1 inclusive (with tolerance of {tolerance}).")
-
- # Check for fill levels within the range [0, 1]
- in_range = np.logical_and(fill_levels > 0, fill_levels < 1)
- in_range_values = fill_levels[in_range]
- in_range_indices = np.where(in_range)
-
- # print("Fill level values within range (0 to 1):")
- # for idx, value in zip(zip(*in_range_indices), in_range_values):
- # if idx == (np.int64(309),):
- # print(f"Index: {idx}, Value: {value}")
-
- # Constants
- pi: float = np.pi # Mathematical constant π
- e: float = 1.60217662e-19 # Elementary charge in Coulombs
- beam_current: float = 3e-3 # Beam current in Amperes
- dwell_time: float = 0.32e-6 # Dwell time in seconds
- A: float = 21.5 # Approximation for Ti6Al4V
- # Ti6Al4V = Ti6Al4V.create_Ti6Al4V(T)
- # A: float = Ti6Al4V.atomic_number
- # print('Ti6Al4V.elements: ', Ti6Al4V.elements)
- # print('Ti6Al4V.composition: ', Ti6Al4V.composition)
- # print('A: ', A)
- dx: float = 3e-6 # Cell size in meters
- dy: float = dx # Cell size in meters
- dz: float = dx # Cell size in meters
- threshold: float = 0 # Beam profile threshold
- # The voltage at detector d can be computed using the transimpedance T (in units of V /A ) of a virtual amplifier in the context of the model
- transimpedance: float = 1e5 # Transimpedance in Volts per Ampere
- eV_threshold = 50 # Threshold energy in eV
- threshold_energy = eV_threshold * e # Threshold energy in joules
- # Define the work distance (height of the electron beam gun)
- # w_d = 12e-3 # Work distance in meters (272 mm)
-
- # Domain setup (Section 2.1.1)
- # Define the 3D domain for the simulation
- x_min, x_max = -0.5 * x_num * dx, 0.5 * x_num * dx # -10 mm to 10 mm
- # print('x_min, x_max:', x_min, x_max)
- y_min, y_max = -0.5 * y_num * dy, 0.5 * y_num * dy # -10 mm to 10 mm
- # print('y_min, y_max:', y_min, y_max)
- z_min, z_max = 0, z_num * dz
- # print('z_min, z_max:', z_min, z_max)
- # x_num = int((x_max - x_min) / dx) + 1
- # y_num = int((y_max - y_min) / dy) + 1
- # z_num = int((z_max - z_min) / dz) + 1
- # print(f"Domain dimensions: x={x_num}, y={y_num}, z={z_num}")
- # domain = np.zeros((z_num, y_num, x_num)) # (z, y , x)
- # print('domain shape: ', np.shape(domain))
- # domain[0:2, :, :] = 1 # Solid cells
- # domain[2, :, :] = 0.5 # Interface layer
- fill_levels_reshaped = np.reshape(fill_levels, (z_num, y_num, x_num))
- # fill_levels_reshaped = fill_levels
- # print('fill_levels_reshaped:', fill_levels_reshaped)
- out_of_range_fill_levels = np.logical_or(fill_levels_reshaped < -tolerance, fill_levels_reshaped > 1 + tolerance)
- if np.any(out_of_range_fill_levels):
- out_of_range_fill_level_values = fill_levels_reshaped[out_of_range_fill_levels]
- out_of_range_fill_level_indices = np.where(out_of_range_fill_levels)
-
- for idx, value in zip(zip(*out_of_range_fill_level_indices), out_of_range_fill_level_values):
- print(f"Index: {idx}, Value: {value}")
- raise ValueError(f"Fill level values must be between 0 and 1 inclusive (with tolerance of {tolerance}).")
-
- sigma = 0.0
- # Different sigma for each dimension
- # sigma = [0.0, 1.0, 0.5] # [z, y, x]
- fill_levels_smoothed = ndimage.gaussian_filter(fill_levels_reshaped, sigma=sigma, mode='nearest', radius=1)
-
- fill_levels_smoothed = np.clip(
- fill_levels_smoothed,
- a_min=np.float64(0.0),
- a_max=np.float64(1.0)
- )
-
- # Check if values are still in bounds
- if not (np.all(fill_levels_smoothed >= 0) and np.all(fill_levels_smoothed <= 1)):
- raise ValueError("Smoothed fill levels out of bounds!")
-
-
- def interpolate_normals(_normal_below, _normal_above, t):
- """Interpolate between two normals using SLERP"""
- # Ensure normals are normalized
- _normal_below = normalize_vector(_normal_below)
- _normal_above = normalize_vector(_normal_above)
-
- # Compute dot product
- dot_product = np.clip(np.dot(_normal_below, _normal_above), -1.0, 1.0)
-
- # If normals are very close, use linear interpolation
- if dot_product > np.cos(np.deg2rad(1)):
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- # If normals are nearly opposite, rotate around arbitrary axis
- elif dot_product < -0.9999:
- # Create perpendicular vector for rotation
- perp = np.array([1.0, 0.0, 0.0]) if abs(_normal_below[1]) < 0.9 else np.array([0.0, 1.0, 0.0])
- perp = normalize_vector(np.cross(_normal_below, perp))
- # Rotate 180 degrees around perp
- _normal_interp = _normal_below * np.cos(np.pi * t) + perp * np.sin(np.pi * t)
- # Otherwise use spherical interpolation
- else:
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + \
- (np.sin(t*theta) / sin_theta) * _normal_above
-
- return normalize_vector(_normal_interp)
-
-
- def interpolate_fill_level_and_normal(_fill_level, surface_normals, _y, _x):
- # Validate input
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
- # Find the indices where the fill level crosses 0.5 in either direction
- indices_above = np.where((_fill_level[:, _y, _x][:-1] > 0.5) & (_fill_level[:, _y, _x][1:] <= 0.5))[0]
- print('indices_above:', indices_above)
- indices_below = np.where((_fill_level[:, _y, _x][:-1] <= 0.5) & (_fill_level[:, _y, _x][1:] > 0.5))[0]
- print('indices_below:', indices_below)
-
- # if len(indices_above) == 0 and len(indices_below) == 0:
- # Fill level does not cross 0.5 at the given (_x, _y) coordinate
- # return None, None, None, None, None, None, None
-
- if len(indices_above) > 0:
- print("len(indices_above) > 0")
- # Fill level crosses 0.5 from above to below
- z_index = indices_above[0]
- print('z_index:', z_index)
- _fill_level_above = _fill_level[z_index, _y, _x]
- print('_fill_level_above:', _fill_level_above)
- _fill_level_below = _fill_level[z_index + 1, _y, _x]
- print('_fill_level_below:', _fill_level_below)
- elif len(indices_below) > 0:
- print("len(indices_below) > 0")
- # Fill level crosses 0.5 from below to above
- z_index = indices_below[0]
- print('z_index:', z_index)
- _fill_level_below = _fill_level[z_index, _y, _x]
- print('_fill_level_below:', _fill_level_below)
- _fill_level_above = _fill_level[z_index + 1, _y, _x]
- print('_fill_level_above:', _fill_level_above)
- else:
- # Fill level does not cross 0.5 at the given (_x, _y) coordinate
- print("No interface cells found")
- return None, None, None, None, None, None, None
-
- # Interpolate the z-coordinate where the fill level is 0.5
- denominator = _fill_level_below - _fill_level_above
- if abs(denominator) < 1e-6:
- # Handle the case where the denominator is very close to zero
- _z_interp = z_index + 0.5
- else:
- _z_interp = z_index + ((0.5 - _fill_level_above) / denominator)
- # _z_interp = z_index + (0.5 - bilinear_interpolate(_x, _y, z_index, _fill_level)) / (bilinear_interpolate(_x, _y, z_index+1, _fill_level) - bilinear_interpolate(_x, _y, z_index, _fill_level))
- print('_z_interp:', _z_interp)
-
- _z_height = (_z_interp + 0.5) * dz
- print('_z_height:', _z_height)
-
- # Interpolate the surface normal at the interpolated z-coordinate
- _normal_below = surface_normals[z_index, _y, _x].astype(np.float64)
- print('_normal_below:', _normal_below)
- _normal_above = surface_normals[z_index + 1, _y, _x].astype(np.float64)
- print('_normal_above:', _normal_above)
-
- t = _z_interp - z_index
-
- # Use spherical interpolation (slerp) for more accurate normal interpolation
- dot_product = np.dot(_normal_below, _normal_above).astype(np.float64)
- if dot_product > np.cos(np.deg2rad(1)): # Threshold for using linear interpolation
- # If normals are very close, use linear interpolation
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- # Use spherical interpolation
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + (np.sin(t*theta) / sin_theta) * _normal_above
- # _normal_interp = quaternion_interpolate(_normal_below, _normal_above, t)
-
- _normal_interp = normalize_vector(_normal_interp)
- print('_normal_interp:', _normal_interp)
-
- return _z_interp, _z_height, _fill_level_below, _fill_level_above, _normal_interp, _normal_below, _normal_above
-
-
- def interpolate_fill_level_and_normal1(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- # Validate input
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at the given (x,y) position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom (reverse order)
- for z in range(len(column)-2, -1, -1): # Start from second-to-last index
- fill_current = column[z]
- fill_next = column[z + 1]
-
- # Check if we cross 0.5 threshold
- if (fill_current - 0.5) * (fill_next - 0.5) <= 0:
- # Found interface
- z_index = z
- _fill_level_above = fill_current
- _fill_level_below = fill_next
-
- print(f"Found interface at z={z_index}")
- print(f"fill_level_above: {_fill_level_above}")
- print(f"fill_level_below: {_fill_level_below}")
-
- # Interpolate z position
- denominator = _fill_level_below - _fill_level_above
- if abs(denominator) < 1e-6:
- _z_interp = z_index + 0.5
- else:
- _z_interp = z_index + ((0.5 - _fill_level_above) / denominator)
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and interpolate normals
- _normal_below = surface_normals[z_index, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z_index + 1, _y, _x].astype(np.float64)
-
- # Interpolate normals
- t = _z_interp - z_index
- dot_product = np.dot(_normal_below, _normal_above)
-
- if dot_product > np.cos(np.deg2rad(1)):
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + \
- (np.sin(t*theta) / sin_theta) * _normal_above
-
- _normal_interp = normalize_vector(_normal_interp)
-
- return _z_interp, _z_height, _fill_level_below, _fill_level_above, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
-
- def interpolate_fill_level_and_normal2(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at the given (x,y) position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom
- for z in range(len(column)-2, -1, -1):
- fill_current = column[z]
- fill_next = column[z + 1]
-
- # Look for significant transition from liquid (>0.5) to gas (<0.2)
- if fill_current > 0.5 and fill_next < 0.2:
- print(f"Found interface at z={z}")
- print(f"fill_level_above: {fill_next}") # gas phase
- print(f"fill_level_below: {fill_current}") # liquid phase
-
- # Interpolate z position
- denominator = fill_next - fill_current
- if abs(denominator) < 1e-6:
- _z_interp = z + 0.5
- else:
- _z_interp = z + ((0.5 - fill_current) / denominator)
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and interpolate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- # Interpolate normal
- t = _z_interp - z
- dot_product = np.dot(_normal_below, _normal_above)
-
- if dot_product > np.cos(np.deg2rad(1)):
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + \
- (np.sin(t*theta) / sin_theta) * _normal_above
-
- _normal_interp = normalize_vector(_normal_interp)
-
- return _z_interp, _z_height, fill_next, fill_current, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
-
- def interpolate_fill_level_and_normal3(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom
- for z in range(len(column)-2, -1, -1):
- current_fill = column[z]
- next_fill = column[z + 1]
-
- # Look for significant liquid-liquid or liquid-gas interface
- if current_fill > 0.5 and next_fill < 0.5:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}") # higher fill level
- print(f"fill_level_above: {next_fill}") # lower fill level
-
- # Interpolate z position
- denominator = next_fill - current_fill
- if abs(denominator) < 1e-6:
- _z_interp = z + 0.5
- else:
- _z_interp = z + ((0.5 - current_fill) / denominator)
-
- print('_z_interp', _z_interp)
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and interpolate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- # Interpolate normal
- t = _z_interp - z
- dot_product = np.dot(_normal_below, _normal_above)
-
- if dot_product > np.cos(np.deg2rad(1)):
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + \
- (np.sin(t*theta) / sin_theta) * _normal_above
-
- _normal_interp = normalize_vector(_normal_interp)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
- def interpolate_fill_level_and_normal4(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom
- for z in range(len(column)-2, -1, -1):
- current_fill = column[z]
- next_fill = column[z + 1]
-
- # Look for transition from solid (=1) to liquid (<1)
- if current_fill == 1.0 and next_fill < 1.0:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {next_fill}") # liquid phase
- print(f"fill_level_above: {current_fill}") # solid phase
-
- # Interpolate z position
- denominator = next_fill - current_fill
- if abs(denominator) < 1e-6:
- _z_interp = z + 0.5
- else:
- _z_interp = z + ((0.5 - current_fill) / denominator)
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and interpolate normals
- _normal_below = surface_normals[z + 1, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z, _y, _x].astype(np.float64)
-
- # Interpolate normal
- t = _z_interp - z
- dot_product = np.dot(_normal_below, _normal_above)
-
- if dot_product > np.cos(np.deg2rad(1)):
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + \
- (np.sin(t*theta) / sin_theta) * _normal_above
-
- _normal_interp = normalize_vector(_normal_interp)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
-
- def interpolate_fill_level_and_normal5(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom
- for z in range(len(column)-2, -1, -1):
- current_fill = column[z]
- next_fill = column[z + 1]
-
- # Look for liquid-gas interface (current > 0, next = 0)
- if current_fill > 0 and next_fill == 0:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}") # liquid phase
- print(f"fill_level_above: {next_fill}") # gas phase
-
- # Interface is at z since this is the last liquid cell
- denominator = next_fill - current_fill
- if abs(denominator) < 1e-6:
- _z_interp = z + 0.5
- else:
- _z_interp = z + ((0.5 - current_fill) / denominator)
-
- print('_z_interp', _z_interp)
- _z_height = (_z_interp + 0.5) * dz
- print("_z_height:", _z_height)
-
- # Get and interpolate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- # Interpolate normal
- t = _z_interp - z
- print("t:", t)
-
- # Use default normal if either normal is invalid
- if np.any(np.isnan(_normal_below)) or np.any(np.isnan(_normal_above)):
- _normal_interp = np.array([0.0, 0.0, 1.0])
- else:
- # Interpolate normal at interface
- _normal_interp = interpolate_normals(_normal_below, _normal_above, t)
-
- # _normal_interp = interpolate_normals(_normal_below, _normal_above, t)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
-
- def compute_surface_normals_sobel(_fill_levels: np.ndarray) -> np.ndarray:
- """
- Compute surface normals using Sobel filter with enhanced numerical precision.
- """
- # Convert input to high precision at the start
- _fill_levels = _fill_levels.astype(np.float64)
-
- # Input validation
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- # Define tolerances for different checks
- TOLERANCE = {
- 'fill_level_bounds': 1e-10,
- 'zero': 1e-10,
- 'norm': 1e-10,
- 'normalization': 1e-8
- }
-
- # Validate domain bounds with high precision
- '''invalid_elements = np.logical_or(_fill_levels < -TOLERANCE['fill_level_bounds'],
- _fill_levels > 1 + TOLERANCE['fill_level_bounds'])
- if np.any(invalid_elements):
- invalid_indices = np.where(invalid_elements)
- for _z, _y, _x in zip(*invalid_indices):
- _value = _fill_levels[_z, _y, _x]
- print(f"Warning: Invalid value {_value} at [z={_z}, y={_y}, x={_x}]")'''
-
- # Compute gradients using Sobel filter with explicit high precision
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64)
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64)
- _grad_z = -1 * ndimage.sobel(_fill_levels, axis=0, output=np.float64) # multiply by -1 to invert the surface normals to point in +z direction
-
- # Stack gradients into normal vector array with maintained precision
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Compute norm with high precision
- norm = np.linalg.norm(_n_S, axis=-1)
-
- # Create masks for valid regions
- interface_mask = (_fill_levels > TOLERANCE['zero']) & (_fill_levels < 1 - TOLERANCE['zero'])
- nonzero_mask = norm > TOLERANCE['norm']
- valid_mask = interface_mask & nonzero_mask
-
- # Normalize vectors only where valid
- # _n_S[valid_mask] /= norm[valid_mask][..., np.newaxis]
- # _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized = np.full_like(_n_S, np.nan)
- _n_S_normalized[valid_mask] = _n_S[valid_mask] / norm[valid_mask, np.newaxis]
-
- # Verify normalization
- if not np.allclose(np.linalg.norm(_n_S_normalized[valid_mask], axis=-1),1.0, atol=TOLERANCE['normalization']):
- raise ValueError("Normalization failed for some vectors")
-
- # return _n_S
- return _n_S_normalized
-
-
- def compute_surface_normals_sobel1(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'fill_level_bounds': 1e-10,
- 'zero': 1e-10,
- 'norm': 1e-10,
- 'normalization': 1e-8
- }
-
- # Compute gradients
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64)
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64)
- _grad_z = -1 * ndimage.sobel(_fill_levels, axis=0, output=np.float64)
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Compute norm
- norm = np.linalg.norm(_n_S, axis=-1)
-
- # Create masks for valid regions:
- # 1. Interface cells (fill level between 0 and 1)
- interface_mask = (_fill_levels > TOLERANCE['zero']) & (_fill_levels < 1 - TOLERANCE['zero'])
-
- # 2. Cells adjacent to fill level transitions
- # filled_mask = _fill_levels >= 1 - TOLERANCE['zero']
- filled_mask = _fill_levels >= 1
- # empty_mask = _fill_levels <= TOLERANCE['zero']
- empty_mask = _fill_levels <= 0
-
- # Check for transitions in z direction
- z_transitions = np.zeros_like(_fill_levels, dtype=bool)
- z_transitions[:-1, :, :] |= (filled_mask[:-1, :, :] & empty_mask[1:, :, :])
- z_transitions[1:, :, :] |= (empty_mask[:-1, :, :] & filled_mask[1:, :, :])
-
- # Combine masks
- valid_mask = interface_mask | z_transitions
-
- # Default normal for sharp transitions (pointing up)
- _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized[..., 2] = 1.0 # Default z-direction normal
-
- # Expand valid_mask to match _n_S shape
- valid_mask_expanded = valid_mask[..., np.newaxis]
-
- # Where we have valid gradients, use them
- gradient_mask = (norm > TOLERANCE['norm'])[..., np.newaxis]
- mask = valid_mask_expanded & gradient_mask
-
- # Normalize vectors where mask is True
- norm_expanded = norm[..., np.newaxis]
- _n_S_normalized = np.where(mask, _n_S / np.where(norm_expanded > 0, norm_expanded, 1.0), _n_S_normalized)
-
- return _n_S_normalized
-
-
- def compute_surface_normals_sobel2(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'fill_level_bounds': 1e-10,
- 'zero': 1e-10,
- 'norm': 1e-10,
- 'normalization': 1e-8
- }
-
- # Compute gradients
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64)
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64)
- _grad_z = -1 * ndimage.sobel(_fill_levels, axis=0, output=np.float64)
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Compute norm
- norm = np.linalg.norm(_n_S, axis=-1)
-
- # Create masks for valid regions:
- # 1. Interface cells (fill level between 0 and 1)
- interface_mask = (_fill_levels > TOLERANCE['zero']) & (_fill_levels < 1 - TOLERANCE['zero'])
-
- # 2. Cells adjacent to fill level transitions
- filled_mask = _fill_levels >= 1 - TOLERANCE['zero']
- empty_mask = _fill_levels <= TOLERANCE['zero']
-
- # Check for transitions in z direction
- z_transitions = np.zeros_like(_fill_levels, dtype=bool)
- z_transitions[:-1, :, :] |= (filled_mask[:-1, :, :] & empty_mask[1:, :, :])
- z_transitions[1:, :, :] |= (empty_mask[:-1, :, :] & filled_mask[1:, :, :])
-
- # Combine masks
- valid_mask = interface_mask | z_transitions
-
- # Initialize output array with NaN
- _n_S_normalized = np.full_like(_n_S, np.nan)
-
- # Create mask for valid gradients
- gradient_mask = norm > TOLERANCE['norm']
-
- # Combine masks and reshape for broadcasting
- mask = valid_mask & gradient_mask
-
- # Normalize only where mask is True
- for i in range(3): # For each component (x, y, z)
- _n_S_normalized[..., i] = np.where(mask,
- _n_S[..., i] / np.where(norm > 0, norm, 1.0),
- np.nan)
-
- return _n_S_normalized
-
-
- def interpolate_normals_new(_normal_below: np.ndarray, _normal_above: np.ndarray, t: float) -> np.ndarray:
- """Interpolate between two normals using SLERP"""
- if not (0 <= t <= 1):
- raise ValueError(f"Interpolation parameter t must be in [0,1], got {t}")
-
- # Ensure normalized vectors
- _normal_below = normalize_vector(_normal_below.astype(np.float64))
- _normal_above = normalize_vector(_normal_above.astype(np.float64))
-
- # Compute dot product
- dot_product = np.clip(np.dot(_normal_below, _normal_above), -1.0, 1.0)
-
- # Thresholds for special cases
- PARALLEL_THRESHOLD = 0.9999
- ANTIPARALLEL_THRESHOLD = -0.9999
-
- # Case 1: Nearly parallel
- if dot_product > PARALLEL_THRESHOLD:
- return normalize_vector((1.0 - t) * _normal_below + t * _normal_above)
-
- # Case 2: Nearly antiparallel
- if dot_product < ANTIPARALLEL_THRESHOLD:
- perp = np.array([1.0, 0.0, 0.0])
- if abs(np.dot(_normal_below, perp)) > 0.9:
- perp = np.array([0.0, 1.0, 0.0])
- perp = normalize_vector(np.cross(_normal_below, perp))
- angle = np.pi * t
- return normalize_vector(_normal_below * np.cos(angle) + perp * np.sin(angle))
-
- # Case 3: Standard SLERP
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
-
- if abs(sin_theta) < 1e-10:
- return normalize_vector(_normal_below)
-
- scale0 = np.sin((1.0 - t) * theta) / sin_theta
- scale1 = np.sin(t * theta) / sin_theta
- return normalize_vector(scale0 * _normal_below + scale1 * _normal_above)
-
-
- def interpolate_fill_level_and_normal_new(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals: {surface_normals.ndim}")
-
- # Get fill levels at position
- _column = _fill_level[:, _y, _x]
-
- # Constants with physical meaning
- INTERFACE_THRESHOLD = 0.5 # Standard VOF interface threshold
- GRADIENT_THRESHOLD = 0.1 # Minimum gradient for interface detection
- TOLERANCE = 1e-6 # Numerical tolerance
-
- # Scan from top to bottom
- for z in range(len(_column)-2, -1, -1):
- current_fill = _column[z]
- next_fill = _column[z + 1]
-
- # Check for significant gradient
- gradient = abs(next_fill - current_fill)
- if gradient > GRADIENT_THRESHOLD:
-
- # Check for interface crossing (VOF method)
- crosses_interface = (
- (current_fill >= (INTERFACE_THRESHOLD - TOLERANCE) and
- next_fill <= (INTERFACE_THRESHOLD + TOLERANCE))
- )
-
- if crosses_interface:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}")
- print(f"fill_level_above: {next_fill}")
-
- denominator = next_fill - current_fill
- if abs(denominator) < TOLERANCE:
- _z_interp = z + 0.5
- t_interface = 0.5
- else:
- # Handle numerical precision near 0.5
- if abs(current_fill - INTERFACE_THRESHOLD) < TOLERANCE:
- t_interface = 0.0
- elif abs(next_fill - INTERFACE_THRESHOLD) < TOLERANCE:
- t_interface = 1.0
- else:
- # Calculate interpolation parameter
- t_interface = (INTERFACE_THRESHOLD - current_fill) / denominator
-
- # Ensure t is in valid range
- if not (0 <= t_interface <= 1):
- error_msg = (
- f"\nInterpolation parameter t_interface={t_interface:.6f} is out of bounds [0,1]\n"
- f"Debug information:\n"
- f" - z-index: {z}\n"
- f" - current_fill: {current_fill:.6f}\n"
- f" - next_fill: {next_fill:.6f}\n"
- f" - denominator: {denominator:.6f}\n"
- f" - INTERFACE_THRESHOLD: {INTERFACE_THRESHOLD}\n"
- f" - Position (y,x): ({_y}, {_x})\n"
- f" - Fill level values around interface:\n"
- )
- # Nearby fill levels for context
- for i in range(max(0, z-2), min(len(_column), z+3)):
- error_msg += f" {i:3d} | {_column[i]:.6f}\n"
-
- raise ValueError(error_msg)
- _z_interp = z + t_interface
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and validate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- _normal_interp = interpolate_normals_new(_normal_below, _normal_above, t_interface)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- return None, None, None, None, None, None, None
-
-
- def compute_surface_normals_sobel_zeros(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with default upward normal
- _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized[..., 2] = 1.0
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
- def compute_surface_normals_sobel_nans(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter with NaN initialization"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with NaN
- _n_S_normalized = np.full_like(_n_S, np.nan)
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
- # Compute surface normals (Section 2.1.1)
- n_S = compute_surface_normals_sobel_nans(fill_levels_smoothed)
- # print('n_S: ', np.shape(n_S))
-
- # Beam profile (Section 2.1)
- # Define the Gaussian beam profile
- sigma: float = 500e-6 / (4 * dx) # Beam profile sigma
- x = np.arange(-np.floor(min(x_num//4, 2*sigma)), 1+np.floor(min((x_num//4), 2*sigma)))
- # print('x:', x)
- y = np.arange(-floor(min(y_num//4, 2*sigma)), 1+floor(min((y_num//4), 2*sigma)))
- # print('y:', y)
- X, Y = np.meshgrid(x, y)
- beam_profile = np.exp(-(X ** 2 + Y ** 2) / (2 * sigma ** 2))
- beam_profile /= np.sum(beam_profile) # Normalize beam profile
- print('beam_profile shape: ', np.shape(beam_profile))
-
- beam_size = np.shape(beam_profile)[0] * np.shape(beam_profile)[1]
- # print('Beam size:', beam_size)
-
- # Detector setup (Section 2.1.1)
- # Define positions of the detectors
- detector_height = 4 * z_max # Detector height in meters (272 mm)
- detector_spacing_x = 1.5 * dx * x_num # Detector spacing in meters (127 mm)
- detector_spacing_y = detector_spacing_x # 1.5 * dy * y_num # Detector spacing in meters (127 mm)
- detector_diameter = 4 * dx # Detector diameter in meters (50 mm)
-
- # Define detector positions relative to the origin
- detector_positions = np.array([
- [0.5 * detector_spacing_x, 0, detector_height],
- [-0.5 * detector_spacing_x, 0, detector_height],
- [0, 0.5 * detector_spacing_y, detector_height],
- [0, -0.5 * detector_spacing_y, detector_height]
- ], dtype=np.float64)
- # print('detector_positions:', detector_positions)
-
-
- # Define detector vertices for solid angle calculation
- def get_detector_vertices(detector_position: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
- # Calculate the side length of the detector
- side_length : float = np.sqrt(3) * (detector_diameter / 2)
- height : float = np.sqrt(3) * (side_length / 2)
-
- # Determine the orientation based on the detector position
- if np.allclose(detector_position[:2], [0, 0]):
- raise ValueError("Detector position cannot be at the origin (0, 0, 0)")
- elif detector_position[0] > 0: # Right detector
- _v1 = detector_position + np.array([2 / 3 * height, 0, 0], dtype=np.float64)
- _v2 = detector_position + np.array([-height / 3, side_length / 2, 0], dtype=np.float64)
- _v3 = detector_position + np.array([-height / 3, -side_length / 2, 0], dtype=np.float64)
- elif detector_position[0] < 0: # Left detector
- _v1 = detector_position + np.array([-2 / 3 * height, 0, 0], dtype=np.float64)
- _v2 = detector_position + np.array([height / 3, -side_length / 2, 0], dtype=np.float64)
- _v3 = detector_position + np.array([height / 3, side_length / 2, 0], dtype=np.float64)
- elif detector_position[1] > 0: # Back detector
- _v1 = detector_position + np.array([0, 2 / 3 * height, 0], dtype=np.float64)
- _v2 = detector_position + np.array([-side_length / 2, -height / 3, 0], dtype=np.float64)
- _v3 = detector_position + np.array([side_length / 2, -height / 3, 0], dtype=np.float64)
- elif detector_position[1] < 0: # Front detector
- _v1 = detector_position + np.array([0, -2 / 3 * height, 0], dtype=np.float64)
- _v2 = detector_position + np.array([side_length / 2, height / 3, 0], dtype=np.float64)
- _v3 = detector_position + np.array([-side_length / 2, height / 3, 0], dtype=np.float64)
- else:
- raise ValueError("Invalid detector position. Must be in one of the four expected positions.")
-
- return _v1, _v2, _v3
-
-
- # Compute solid angle using Equation A1 (Appendix A)
- def compute_solid_angle(vector1: np.ndarray, vector2: np.ndarray, vector3: np.ndarray) -> float:
- # Input validation
- if not isinstance(vector1, np.ndarray) or not isinstance(vector2, np.ndarray) or not isinstance(vector3, np.ndarray):
- raise ValueError("Inputs must be numpy arrays")
- if vector1.shape != (3,) or vector2.shape != (3,) or vector3.shape != (3,):
- raise ValueError("Inputs must be 3D vectors")
- if not np.isfinite(vector1).all() or not np.isfinite(vector2).all() or not np.isfinite(vector3).all():
- raise ValueError("Inputs must be finite")
-
- # Convert to high precision once
- vectors = [vec.astype(np.float64) for vec in (vector1, vector2, vector3)]
- cross_product = np.cross(vectors[1], vectors[2])
- # Check if the cross product is nearly zero
- if np.linalg.norm(cross_product) < 1e-16:
- raise ValueError("Vectors are nearly coplanar")
-
- # Compute solid angle
- try:
- # Compute the numerator of the solid angle formula
- numerator: float = np.dot(vectors[0], cross_product)
- # Check if the numerator is nearly zero
- if np.abs(numerator) < 1e-16:
- raise ValueError("Vectors are nearly parallel")
-
- norms = [np.linalg.norm(v) for v in vectors]
-
- # Compute the denominator of the solid angle formula
- denominator = (
- norms[0] * norms[1] * norms[2] +
- np.dot(vectors[0], vectors[1]) * norms[2] +
- np.dot(vectors[0], vectors[2]) * norms[1] +
- np.dot(vectors[1], vectors[2]) * norms[0]
- )
-
- # Compute the solid angle in radians
- solid_angle = 2 * np.arctan2(numerator, denominator)
- # print("Solid angle:", solid_angle)
- return solid_angle
- except ZeroDivisionError:
- raise ValueError("Error computing solid angle: division by zero")
- except np.linalg.LinAlgError:
- raise ValueError("Error computing solid angle: linear algebra error")
-
-
- # BSE coefficient functions (Equations 11 and 10, Section 2.1.2)
- def eta_0(_A: float) -> float:
- """Equation 11 from the paper"""
- _eta_0: float = -0.0254 + 0.016 * _A - 1.86e-4 * _A ** 2 + 8.3e-7 * _A ** 3
- print('eta_0: ', _eta_0)
- return _eta_0
-
-
- # BSE Coefficient for single atomic targets, found to hold in the range of 10–100 keV
- def eta(theta: float, _A: float) -> float:
- """Equation 10 from the paper"""
- B: float = 0.89
- _eta: float = B * (eta_0(_A) / B) ** np.cos(theta)
- print('eta: ', _eta)
- return _eta
-
-
- '''
- # Compute the weighted average BSE coefficient for the alloy
- def compute_weighted_eta(theta: float, alloy: Alloy) -> float:
- """Equation 12 from the paper"""
- weighted_eta: float = 0
- for element, weight_fraction in zip(alloy.elements, alloy.composition):
- # print('element: ', element)
- # print('weight_fraction: ', weight_fraction)
- _A: float = element.atomic_number
- # print('A: ', _A)
- eta_i: float = eta(theta, _A)
- # print('eta_i: ', eta_i)
- weighted_eta += weight_fraction * eta_i
- # print('weighted_eta: ', weighted_eta)
- return weighted_eta
- '''
-
-
- # Correction function C(theta) (Equation 14, Section 2.1.3)
- def C(theta: float) -> float:
- print("theta1:", theta)
- """Equation 14 from the paper"""
- theta = np.rad2deg(theta) # Ensure angle is in degrees for C function
- _C: float = 1.01 - 4.06 * theta + 1.36e-5 * theta ** 2 - 1.71e-6 * theta ** 3
- print('C: ', _C)
- # Check if C is negative and raise an exception
- if _C < 0:
- raise ValueError(f"C function returned a negative value for theta: {theta}. "
- f"Please check the input and the C function implementation.")
- return _C
-
-
- # Scattering function g_BSE (Equation 13, Section 2.1.3)
- def g_BSE(_theta: float, _alpha: float, _beta: float) -> float:
- print("_theta, _alpha, _beta:", _theta, _alpha, _beta)
- """Equation 13 from the paper"""
- # k: float = np.rad2deg(_theta) / 13 # Correct implementation as per the paper
- # print("k: ", k)
- diffusive_part: float = (eta_0(A) / pi) * np.cos(_alpha)
- print(np.cos(_alpha))
- print("diffusive_part:", diffusive_part)
- # reflective_part: float = ((eta(_theta, A) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
- # if np.abs(reflective_part) < 1e-10:
- # reflective_part = 0
- # print("reflective_part:", reflective_part)
-
- # Compute a weighted average BSE coefficient for the alloy
- # reflective_part: float = ((compute_weighted_eta(_theta, Ti6Al4V) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
-
- # print('g_BSE: ', diffusive_part + reflective_part)
- return diffusive_part #+ reflective_part
-
-
- # Calculate the beam direction vector p_E
- def calculate_p_E(_beam_x, _beam_y):
- print(f"calculate_p_E({_beam_x}, {_beam_y}):")
- # Calculate the vector from the beam gun to the beam location
- beam_vector = np.array([_beam_x * dx, _beam_y * dy, -detector_height], dtype=np.float64)
- return normalize_vector(beam_vector)
-
- def calculate_p_E1(_beam_x, _beam_y, _height):
- """Calculate normalized beam direction vector from gun to surface point
-
- Args:
- _beam_x: Beam x-position relative to center (in cells)
- _beam_y: Beam y-position relative to center (in cells)
- _height: Height difference between gun and surface point (m)
- """
- print(f"calculate_p_E({_beam_x}, {_beam_y}):")
- print("_height:", _height)
-
- # Calculate vector components
- x_component = _beam_x * dx # Convert cells to meters
- y_component = _beam_y * dy # Convert cells to meters
- z_component = _height # Already in meters
-
- # Create beam vector (from gun to surface point)
- beam_vector = np.array([x_component, y_component, z_component], dtype=np.float64)
-
- # Verify downward direction
- if normalize_vector(beam_vector)[2] > 0:
- raise ValueError("Beam vector pointing upward")
-
- # Normalize and return
- return normalize_vector(beam_vector)
-
-
- # Calculate phi as the angle between p_E and the negative z-axis
- def calculate_phi(_p_E) -> float:
- z_axis = np.array([0.0, 0.0, -1.0], dtype=np.float64)
-
- # Check if _p_E is approximately equal to the negative z-axis
- '''if np.allclose(_p_E, z_axis, rtol=1e-8, atol=1e-8):
- _phi_rad = 0.0
- else:
- # Calculate the angle between _p_E and the z-axis using arccos
- _phi_rad = calculate_angle(_p_E, z_axis)'''
- _phi_rad = calculate_angle(_p_E, z_axis)
- print(f"calculate_phi({_p_E}, {_phi_rad}):")
- return _phi_rad
-
-
- def set_beam_movement(_x_center, _y_center, _x_offset, _y_offset, _movement_type):
- if _movement_type == 'X':
- if _x_offset >= 0:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset - 1
- else:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset + 1
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center + _y_offset + 1
- elif _movement_type == 'Y':
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center + _x_offset + 1
- if _y_offset >= 0:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset - 1
- else:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset + 1
- elif _movement_type == 'XY':
- if _x_offset >= 0:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset - 1
- else:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset + 1
- if _y_offset >= 0:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset - 1
- else:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset + 1
- else:
- raise ValueError("movement_type must be 'X' or 'Y' or 'XY'")
-
- return _beam_x_start, _beam_x_end, _beam_y_start, _beam_y_end
-
-
- # Initialize arrays to store electron counts
- electrons_per_detector_per_step = np.zeros((x_num, y_num, len(detector_positions)))
- total_electrons_per_step = np.zeros((x_num, y_num))
-
- # Additional arrays to capture detailed information
- electrons_per_detector_per_beam_profile = np.zeros(
- (beam_profile.shape[0], beam_profile.shape[1], len(detector_positions)))
- total_electrons_per_beam_profile = np.zeros((beam_profile.shape[0], beam_profile.shape[1]))
-
- # Calculate the center indices
- x_center = x_num // 2
- y_center = y_num // 2
- # print('x_center, y_center:', x_center, y_center)
-
- x_offset: int = 0
- y_offset: int = 0
- movement_type: str = 'Y'
-
- # Beam X-Movement
- # beam_x_start, beam_x_end = x_center - x_offset, x_center + x_offset + 1
- # beam_y_start, beam_y_end = y_center + y_offset, y_center + y_offset + 1
-
- # Beam Y-Movement
- # beam_x_start, beam_x_end = x_center + x_offset, x_center + x_offset + 1
- # beam_y_start, beam_y_end = y_center - y_offset, y_center + y_offset + 1
-
- # For X or Y Movement
- beam_x_start, beam_x_end, beam_y_start, beam_y_end = set_beam_movement(
- x_center, y_center, x_offset, y_offset, movement_type)
-
- # Determine the step size based on the direction of movement
- step_x = 1 if beam_x_start <= beam_x_end else -1
- step_y = 1 if beam_y_start <= beam_y_end else -1
-
- # print(beam_x_start, beam_x_end, beam_y_start, beam_y_end, step_x, step_y)
-
-
- z_heights = np.zeros(beam_size)
- fill_level_b = np.zeros(beam_size)
- fill_level_a = np.zeros(beam_size)
- normals = np.zeros((beam_size, 3))
- normals_b = np.zeros((beam_size, 3))
- normals_a = np.zeros((beam_size, 3))
- cell_counter = 0
-
-
- def normalize_vector(vector: np.ndarray) -> np.ndarray:
- """Normalize a 3D vector."""
- vector = vector.astype(np.float64) # Convert to float64
- norm = np.linalg.norm(vector)
- if norm < 1e-18: # Avoid division by zero
- return vector # Return original if nearly zero
- return vector / norm # Normalize
-
-
- def calculate_angle(_v1: np.ndarray, _v2: np.ndarray) -> float:
- """Calculate the angle between two normalized vectors."""
- _v1 = normalize_vector(_v1)
- _v2 = normalize_vector(_v2)
-
- dot_product = np.clip(np.dot(_v1, _v2), -1.0, 1.0)
-
- return np.arccos(dot_product) # Return angle in radians
-
- # Main simulation loop to iterate over the simulation domain
- for beam_x in range(beam_x_start, beam_x_end, step_x):
- for beam_y in range(beam_y_start, beam_y_end, step_y):
- beam_loc = (beam_x, beam_y)
- print(f"Beam center is at domain location: {beam_loc}")
-
- # Calculate the beam direction vector p_E
- # >> p_E = calculate_p_E(beam_x - x_center, beam_y - y_center)
- # >> print('p_E: ', p_E)
-
- # Calculate phi as the angle between p_E and the negative z-axis
- # phi_rad = np.arctan2(np.sqrt(dx**2 + dy**2), w_d)
- # >> phi_rad = calculate_phi(p_E)
- # print(f"phi_rad: {phi_rad}")
-
- # Reset counters for this beam location
- total_electrons: int = 0
- electrons_per_detector = np.zeros(len(detector_positions))
-
- # Iterate over the beam profile
- for profile_x in range(beam_profile.shape[0]):
- # for profile_x in [beam_profile.shape[0]//2 + 0]:
- for profile_y in range(beam_profile.shape[1]):
- # for profile_y in [beam_profile.shape[1]//2]:
- print(f"Beam profile shape ({profile_x}, {profile_y})")
- domain_x = beam_loc[0] - beam_profile.shape[0] // 2 + profile_x
- domain_y = beam_loc[1] - beam_profile.shape[1] // 2 + profile_y
- print(f"Domain ({domain_x}, {domain_y})")
-
- # p_E = calculate_p_E(domain_x - x_center, domain_y - y_center)
- # print('p_E: ', p_E)
- # phi_rad = calculate_phi(p_E)
-
- if (0 <= domain_x < x_num) and (0 <= domain_y < y_num):
- if beam_profile[profile_x, profile_y] > threshold:
- print('finding interface_z_indices')
- # print('interface_z_indices:', interface_z)
- z_interp, z_height, fill_level_below, fill_level_above, normal_interp, normals_below, normals_above \
- = interpolate_fill_level_and_normal_new(fill_levels_smoothed, n_S, domain_y, domain_x)
- if z_interp is not None:
- z_heights[cell_counter] = z_height
- fill_level_b[cell_counter] = fill_level_below
- fill_level_a[cell_counter] = fill_level_above
- normals[cell_counter] = normal_interp
- normals_b[cell_counter] = normals_below
- normals_a[cell_counter] = normals_above
- cell_counter += 1
- n_S_local = normal_interp
- print('n_S_local:', n_S_local)
- if not np.allclose(np.linalg.norm(n_S_local), 1.0, atol=1e-5):
- raise ValueError(f"n_S_local is not normalized, {np.linalg.norm(n_S_local)}")
-
- p_E = calculate_p_E1(domain_x - x_center, domain_y - y_center, -(detector_height-z_height))
- print('p_E: ', p_E)
-
- theta_local = calculate_angle(-p_E, n_S_local)
- print('theta_local:', theta_local)
- # if theta_local < 0.0 or theta_local > 0.011:
- # raise ValueError(f"theta_local value should be 0.0, found {theta_local}")
- # theta_local: 0.016
-
- scattering_point = np.array([
- (domain_x - beam_x) * dx,
- (domain_y - beam_y) * dy,
- z_height
- ], dtype=np.float64)
- print(f'scattering_point ({domain_x}, {domain_y}, {z_interp}):', scattering_point)
-
- for det_idx, det_pos in enumerate(detector_positions):
- # print(f"Detector {det_idx + 1} position: {det_pos}")
- det_idx: int
- # det_pos: np.ndarray[float]
- print(f'det_pos[{det_idx + 1}]: {det_pos}')
-
- # Calculate the vector from the scattering point to the detector
- det_vec = det_pos - scattering_point
- det_vec = normalize_vector(det_vec)
- print(f'det_vec[{det_idx + 1}]: {det_vec}')
- if not np.allclose(np.linalg.norm(det_vec), 1.0):
- raise ValueError("det_vec is not normalized.")
-
- # alpha = np.arccos(np.clip(np.dot(det_vec, n_S_local), -1.0, 1.0))
- alpha = calculate_angle(det_vec, n_S_local)
- print(f'alpha: {alpha}')
-
- if alpha > np.pi / 2:
- continue
-
- p_E_reflected = p_E - 2.0 * np.dot(p_E, n_S_local) * n_S_local
- # p_E_reflected = np.where(np.abs(p_E_reflected) < 1e-8, 0.0, p_E_reflected)
- p_E_reflected = normalize_vector(p_E_reflected)
- print(f'p_E_reflected: {p_E_reflected}')
-
- # beta = np.arccos(np.clip(np.dot(p_E_reflected, det_vec), -1.0, 1.0))
- beta = calculate_angle(p_E_reflected, det_vec)
- print(f'beta: {beta}')
-
- # Apply energy threshold
- # electron_energy = N_0 * e * g_bse_value
- # print('electron_energy: ', electron_energy)
-
- # Calculate the initial energy of the primary electrons (E_0)
- V_acc = 60e3 # Accelerating voltage in Volts
- E_0 = V_acc * e # Initial energy of primary electrons in joules
-
- # Calculate the energy of the backscatter electrons
- E_BSE = E_0 * eta(theta_local, A)
- # Compute a weighted average BSE coefficient for the alloy
- # E_BSE = E_0 * compute_weighted_eta(theta_local, Ti6Al4V)
- # print('E_BSE: ', E_BSE)
-
- # Bifurcation logic for PE and SE
- # if E_BSE > threshold_energy:
- # Backscatter electrons (BSE)
- # Compute scattering function g_bse_value
- g_bse_value = g_BSE(theta_local, alpha, beta)
- print(f'g_bse_value: {g_bse_value}')
-
- # Ensure g_bse_value is non-negative
- if g_bse_value <= 0:
- raise ValueError(f"g_bse_value is {g_bse_value} (zero negative), check calculations.")
-
- # Compute solid angle (Equation A1)
- v1, v2, v3 = get_detector_vertices(det_pos)
- normal = np.cross(v2 - v1, v3 - v1)
- normal = normalize_vector(normal)
-
- dOmega = compute_solid_angle(v1 - scattering_point,
- v2 - scattering_point,
- v3 - scattering_point)
- print(f"Detector {det_idx + 1} solid angle: {dOmega}")
-
- # Compute number of electrons collected by the detector (Equation A2)
- # The beam ray stays on each grid point ij for the dwell time t_d and distributes N_0 PE into it
- N_0 = beam_current * dwell_time / e # Equation 3, Section 2.1.1
- print("N_0:", N_0)
-
- # There are N1 electrons in the chamber after the 1st scattering event
- N_1 = N_0 * eta(theta_local, A) # Equation 4, Section 2.1.1
- # Compute a weighted average BSE coefficient for the alloy
- # N_1 = N_0 * compute_weighted_eta(theta_local, Ti6Al4V) # Equation 4, Section 2.1.1
- print("N_1:", N_1)
-
- # The number of electrons N_1_d collected by a general detector d
- N_1_d = N_1 * np.sum(g_bse_value * dOmega) # Equation 5, Section 2.1.1
-
- # Accumulate electrons per detector for this beam profile point
- electrons_per_detector_per_beam_profile[profile_x, profile_y, det_idx] += N_1_d
-
- # Accumulate electrons per detector for this domain point
- electrons_per_detector_per_step[beam_x, beam_y, det_idx] += N_1_d
-
- # Add to total electrons for this beam profile point
- total_electrons_per_beam_profile[profile_x, profile_y] += N_1_d
-
- # Add to total electrons for this domain point
- total_electrons_per_step[beam_x, beam_y] += N_1_d
-
- # Accumulate electrons per detector
- electrons_per_detector[det_idx] += N_1_d
-
- # Add to total electrons
- total_electrons += N_1_d
-
- # print(f"After addition: electrons_per_detector[{det_idx + 1}]: {electrons_per_detector[det_idx]}")
- else:
- print(f"No interface found at ({domain_x}, {domain_y})")
-
- # Print detailed information for this beam profile point
- # print(f"Electrons per detector at this beam profile point ({profile_x}, {profile_y}): {electrons_per_detector_per_beam_profile[profile_x, profile_y]}")
- # print(f"Total electrons at this beam profile point ({profile_x}, {profile_y}): {total_electrons_per_beam_profile[profile_x, profile_y]}")
-
- # Print detailed information for this domain point after iterating over the beam profile
- # print(f"Electrons per detector at this domain point ({beam_x}, {beam_y}): {electrons_per_detector_per_step[beam_x, beam_y]}")
- # print(f"Total electrons at this domain point ({beam_x}, {beam_y}): {total_electrons_per_step[beam_x, beam_y]}")
-
- # print("Shape of electrons_per_detector_per_step:", electrons_per_detector_per_step.shape)
- # print("Number of lines:", len(lines))
-
- # Compute final results
- total_electrons_all_steps = np.sum(total_electrons_per_step)
- total_electrons_per_detector = np.sum(electrons_per_detector_per_step, axis=(0, 1))
-
- # Print final results
- # print(f"Total electrons hitting all detectors: {total_electrons_all_steps}")
- # for i, electrons in enumerate(total_electrons_per_detector):
- # print(f"Total electrons hitting detector {i + 1}: {electrons}")
-
- # for beam_x in range(beam_x_start, beam_x_end, step_x):
- # for beam_y in range(beam_y_start, beam_y_end, step_y):
- # print(f"Electrons per detector at this domain point ({beam_x}, {beam_y}): {electrons_per_detector_per_step[beam_x, beam_y]}")
- # print('end simulation')
- # print(electrons_per_detector_per_step)
-
- # print("total_electrons:", total_electrons, "at timestep", current_timestep)
- # print("electrons_per_det:", electrons_per_det, "at timestep", current_timestep)
- # print("electrons_per_beam_loc:", electrons_per_beam_loc, "at timestep", current_timestep)
- # print("electrons_per_beam_loc_per_det:", electrons_per_beam_loc_per_det, "at timestep", current_timestep)
- return electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, normals, normals_b, normals_a
-
-def read_fill_levels(filepath):
- """Read only the FillingFrac grid from VDB file"""
- try:
- # Try to read just the FillingFrac grid
- grid = vdb.read(filepath, "FillingFrac")
-
- # Get dimensions
- dims = grid.evalActiveVoxelDim()
- print(f"Grid dimensions: {dims}")
-
- # Create array
- fill_array = np.zeros(dims, dtype=np.float64)
-
- # Copy data
- grid.copyToArray(fill_array)
-
- # Transpose to ray tracer's expected order (z, y, x)
- fill_array = np.transpose(fill_array, (2, 1, 0))
-
- # Create standardized array with 256 z cells
- standardized_array = np.zeros((256, dims[1], dims[0]), dtype=np.float64)
-
- # Copy data from original array to standardized array
- standardized_array[:fill_array.shape[0], :, :] = fill_array
-
- # Apply Gaussian smoothing
- sigma = 1.0
- standardized_array_smoothed = ndimage.gaussian_filter(standardized_array, sigma=sigma, mode='nearest', radius=1)
-
- # Print values
- '''x = 157
- y = 196
- print(f"\nFill levels at ({x}, {y}):")
- print("z_index | fill_level")
- print("-" * 20)
-
- # Get fill levels
- fill_levels = standardized_array_smoothed[:, y, x]
- valid_levels = fill_levels[(fill_levels > 0) & (fill_levels < 1)]
-
- if len(valid_levels) > 0:
- for z, val in enumerate(fill_levels):
- if 0 <= val <= 1:
- print(f"{z:7d} | {val:.6f}")
- else:
- print("No fill level found between 0 and 1")'''
-
- return standardized_array_smoothed
-
- except Exception as e:
- print(f"Error reading {filepath}")
- print(f"Error details: {str(e)}")
- try:
- metadata = vdb.readMetadata(filepath)
- print("\nFile metadata:")
- for key, value in metadata.items():
- print(f"{key}: {value}")
- except Exception as meta_e:
- print(f"Error reading metadata: {str(meta_e)}")
- return None
-
-def process_single_timestep(filepath, timestep):
- """Process a single VDB file and run ray tracer"""
- try:
- # Read the FillingFrac grid
- grid = vdb.read(filepath, "FillingFrac")
-
- # Get dimensions
- dims = grid.evalActiveVoxelDim()
- x_num, y_num = dims[0], dims[1]
- z_num = 256 # dims[2] # Fixed z dimension
- print(f"Domain dimensions: {x_num} x {y_num} x {z_num}")
-
- # Read fill levels
- fill_levels = read_fill_levels(filepath)
- if fill_levels is None:
- raise ValueError("Could not read fill levels from VDB file")
-
- '''if fill_levels is not None:
- print("\nSuccessfully read fill levels")
- print(f"Array shape: {np.shape(fill_levels)}")
-
- try:
- # Convert to numpy array if not already
- fill_array = np.asarray(fill_levels)
-
- # Get non-zero elements
- non_zero = fill_array[fill_array > 0]
-
- # Get interface cells (0 < fill_level < 1)
- interface_cells = fill_array[(fill_array > 0) & (fill_array < 1)]
-
- # Print statistics
- if len(non_zero) > 0:
- print("\nNon-zero cells statistics:")
- print(f"Number of non-zero cells: {len(non_zero)}")
- print(f"Min non-zero value: {np.min(non_zero):.6f}")
- print(f"Max value: {np.max(fill_array):.6f}")
-
- if len(interface_cells) > 0:
- print("\nInterface cells statistics:")
- interface_percentage = (len(interface_cells) * 100.0) / len(non_zero)
- print(f"Number of interface cells: {len(interface_cells)} ({interface_percentage:.3f}%)")
- print(f"Min interface value: {np.min(interface_cells):.6f}")
- print(f"Max interface value: {np.max(interface_cells):.6f}")
- print(f"Mean interface value: {np.mean(interface_cells):.6f}")
- else:
- print("\nNo interface cells found (no fill levels between 0 and 1)")
-
- except Exception as e:
- print(f"Error processing array: {str(e)}")'''
-
- # Run ray tracer
- print(f"\nProcessing timestep {timestep}")
- electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, \
- normals, normals_b, normals_a = runRayTracer(
- x_num,
- y_num,
- z_num,
- # fill_levels.flatten(),
- fill_levels,
- timestep
- )
-
- return electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, \
- normals, normals_b, normals_a
-
- except Exception as e:
- print(f"Error processing timestep {timestep}: {str(e)}")
- raise
-
-def main():
- """Main function to process VDB files and run ray tracer"""
- # Directory containing VDB files
- directory = "/local/ca00xebo/softwares/KiSSAM/Markl-plate"
-
- # Get sorted list of VDB files
- vdb_files = []
- for f in os.listdir(directory):
- if f.startswith('geometry-0000') and f.endswith('.vdb'):
- try:
- # Extract timestep number from filename
- timestep = int(f.split('-')[1].split('.')[0])
- vdb_files.append((timestep, f))
- except (IndexError, ValueError) as e:
- print(f"Skipping file {f}: {str(e)}")
-
- # Sort by timestep
- vdb_files.sort() # sort based on timestep number
-
- # Slice the first 120 files
- vdb_files = vdb_files[:120]
-
- total_timesteps = len(vdb_files)
- print(f"Found {total_timesteps} VDB files to process")
- electrons_across_timesteps = np.zeros((total_timesteps, 4))
-
- # Process each timestep
- _results = []
- for timestep, vdb_file in enumerate(vdb_files):
- filepath = os.path.join(directory, vdb_file[1]) # vdb_file[1] contains the filename
- try:
- print(f"\nProcessing file {vdb_file[1]} (timestep {timestep})")
-
- # Process this timestep
- _result = process_single_timestep(filepath, timestep)
-
- # Debug electron counts
- '''print("\nElectron count debug:")
- print(f"Shape of electrons_per_detector_per_step: {_result[0].shape}")
- print(f"Min value: {np.min(_result[0])}")
- print(f"Max value: {np.max(_result[0])}")
- print(f"Mean value: {np.mean(_result[0])}")'''
-
- # Store electron counts
- total_electrons_for_detectors = np.sum(_result[0], axis=(0, 1))
- print(f"Total electrons before assignment: {total_electrons_for_detectors}")
-
- electrons_across_timesteps[timestep] = total_electrons_for_detectors
- print(f"Electrons for timestep {timestep}:")
- print(f"Raw values: {electrons_across_timesteps[timestep]}")
- # print(f"Sum: {np.sum(electrons_across_timesteps[timestep])}")
- print("electrons_across_timesteps:", electrons_across_timesteps)
-
- _results.append(_result)
- print(f"Completed timestep {timestep}")
-
- except Exception as e:
- print(f"Error processing timestep {timestep}: {str(e)}")
- raise
-
- # Final debug output
- '''print("\nFinal electrons_across_timesteps:")
- print(f"Shape: {electrons_across_timesteps.shape}")
- print(f"Min value: {np.min(electrons_across_timesteps)}")
- print(f"Max value: {np.max(electrons_across_timesteps)}")
- print(f"Mean value: {np.mean(electrons_across_timesteps)}")'''
-
- # Plot results
- plot_detector_results1(electrons_across_timesteps)
-
- return _results
-
-def plot_detector_results(electrons_per_detector_time):
- """Plot electron counts for each detector across timesteps"""
- plt.figure(figsize=(12, 8))
- timesteps = range(len(electrons_per_detector_time))
-
- # Plot for each detector
- for detector in range(4):
- plt.plot(timesteps, electrons_per_detector_time[:, detector],
- label=f'Detector {detector+1}',
- marker='o')
-
- plt.title("Electrons Captured by Detectors Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Number of Electrons")
- plt.legend()
- plt.grid(True)
-
- # Add value annotations
- for detector in range(4):
- for t, value in zip(timesteps, electrons_per_detector_time[:, detector]):
- plt.annotate(f'{value:.2e}',
- (t, value),
- textcoords="offset points",
- xytext=(0,10),
- ha='center',
- rotation=45)
-
- plt.tight_layout()
- plt.savefig("detector_electrons_over_time(VDB-120).png", dpi=300, bbox_inches='tight')
- plt.close()
-
- # Create additional visualization - normalized values
- plt.figure(figsize=(12, 8))
- normalized_electrons = electrons_per_detector_time / np.max(electrons_per_detector_time)
-
- for detector in range(4):
- plt.plot(timesteps, normalized_electrons[:, detector],
- label=f'Detector {detector+1}',
- marker='o')
-
- plt.title("Normalized Electrons Captured by Detectors Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Normalized Electron Count")
- plt.legend()
- plt.grid(True)
- plt.tight_layout()
- plt.savefig("detector_electrons_over_time_normalized(VDB-120).png", dpi=300, bbox_inches='tight')
- plt.close()
-
-def plot_detector_results1(electrons_per_detector_time):
- """
- Plot electron counts for detector pairs across timesteps
-
- Args:
- electrons_per_detector_time: Array of shape (num_timesteps, 4) containing
- electron counts for each detector at each timestep
- :param electrons_per_detector_time:
- """
- # Create figure for detectors 1 and 2
- plt.figure(figsize=(12, 8))
- timesteps = range(len(electrons_per_detector_time))
-
- # Plot detectors 1 and 2
- plt.plot(timesteps, electrons_per_detector_time[:, 0],
- label='Detector 1', marker='o', linewidth=1, markevery=20)
- plt.plot(timesteps, electrons_per_detector_time[:, 1],
- label='Detector 2', marker='s', linewidth=1, markevery=20)
-
- # Add value annotations
- '''for t in timesteps:
- plt.annotate(f'{electrons_per_detector_time[t, 0]:.2e}',
- (t, electrons_per_detector_time[t, 0]),
- textcoords="offset points",
- xytext=(0,10),
- ha='center',
- rotation=45)
- plt.annotate(f'{electrons_per_detector_time[t, 1]:.2e}',
- (t, electrons_per_detector_time[t, 1]),
- textcoords="offset points",
- xytext=(0,-15),
- ha='center',
- rotation=45)'''
-
- plt.title("Electrons Captured by Detectors 1 and 2 Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Number of Electrons")
- plt.legend()
- plt.grid(True)
- plt.tight_layout()
- plt.savefig("detectors_1_2_over_time(VDB-120).png", dpi=300, bbox_inches='tight')
- plt.close()
-
- # Create figure for detectors 3 and 4
- plt.figure(figsize=(12, 8))
-
- # Plot detectors 3 and 4
- plt.plot(timesteps, electrons_per_detector_time[:, 2],
- label='Detector 3', marker='o', linewidth=1, markevery=20)
- plt.plot(timesteps, electrons_per_detector_time[:, 3],
- label='Detector 4', marker='s', linewidth=1, markevery=20)
-
- # Add value annotations
- '''for t in timesteps:
- plt.annotate(f'{electrons_per_detector_time[t, 2]:.2e}',
- (t, electrons_per_detector_time[t, 2]),
- textcoords="offset points",
- xytext=(0,10),
- ha='center',
- rotation=45)
- plt.annotate(f'{electrons_per_detector_time[t, 3]:.2e}',
- (t, electrons_per_detector_time[t, 3]),
- textcoords="offset points",
- xytext=(0,-15),
- ha='center',
- rotation=45)'''
-
- plt.title("Electrons Captured by Detectors 3 and 4 Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Number of Electrons")
- plt.legend()
- plt.grid(True)
- plt.tight_layout()
- plt.savefig("detectors_3_4_over_time(VDB-120).png", dpi=300, bbox_inches='tight')
- plt.close()
-
-if __name__ == "__main__":
- results = main()
\ No newline at end of file
diff --git a/apps/showcases/FreeSurface/raytracerVDB-multi-000003.py b/apps/showcases/FreeSurface/raytracerVDB-multi-000003.py
deleted file mode 100644
index 5902838d7457a19307169eda7736849e41a0c48a..0000000000000000000000000000000000000000
--- a/apps/showcases/FreeSurface/raytracerVDB-multi-000003.py
+++ /dev/null
@@ -1,1797 +0,0 @@
-import pyopenvdb as vdb
-import os
-from math import floor
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib.colors as colors
-from matplotlib.animation import FuncAnimation
-from typing import Tuple, List
-import scipy.ndimage as ndimage
-from scipy.ndimage import gaussian_filter
-from scipy.spatial.transform import Rotation as Ro
-from scipy.spatial.transform import Slerp
-from scipy.interpolate import PchipInterpolator
-from sympy.printing.tree import print_node
-
-
-def runRayTracer(x_num, y_num, z_num, fill_levels, current_timestep):
-
- print("current_timestep", current_timestep)
- # fill_levels = np.clip(fill_levels, 0, 1)
- fill_levels = fill_levels.astype(np.float64)
-
- # Define the tolerance
- tolerance = 1e-8
-
- # Check if any value in fill_levels is outside the [0, 1] range
- out_of_range = np.logical_or(fill_levels < -tolerance, fill_levels > 1 + tolerance)
- if np.any(out_of_range):
- out_of_range_values = fill_levels[out_of_range]
- out_of_range_indices = np.where(out_of_range)
-
- for idx, value in zip(zip(*out_of_range_indices), out_of_range_values):
- print(f"Index: {idx}, Value: {value}")
-
- # raise ValueError(f"Fill level values must be between 0 and 1 inclusive (with tolerance of {tolerance}).")
-
- # Check for fill levels within the range [0, 1]
- in_range = np.logical_and(fill_levels > 0, fill_levels < 1)
- in_range_values = fill_levels[in_range]
- in_range_indices = np.where(in_range)
-
- # print("Fill level values within range (0 to 1):")
- # for idx, value in zip(zip(*in_range_indices), in_range_values):
- # if idx == (np.int64(309),):
- # print(f"Index: {idx}, Value: {value}")
-
- # Constants
- pi: float = np.pi # Mathematical constant π
- e: float = 1.60217662e-19 # Elementary charge in Coulombs
- beam_current: float = 3e-3 # Beam current in Amperes
- dwell_time: float = 0.32e-6 # Dwell time in seconds
- A: float = 21.5 # Approximation for Ti6Al4V
- # Ti6Al4V = Ti6Al4V.create_Ti6Al4V(T)
- # A: float = Ti6Al4V.atomic_number
- # print('Ti6Al4V.elements: ', Ti6Al4V.elements)
- # print('Ti6Al4V.composition: ', Ti6Al4V.composition)
- # print('A: ', A)
- dx: float = 3e-6 # Cell size in meters
- dy: float = dx # Cell size in meters
- dz: float = dx # Cell size in meters
- threshold: float = 0 # Beam profile threshold
- # The voltage at detector d can be computed using the transimpedance T (in units of V /A ) of a virtual amplifier in the context of the model
- transimpedance: float = 1e5 # Transimpedance in Volts per Ampere
- eV_threshold = 50 # Threshold energy in eV
- threshold_energy = eV_threshold * e # Threshold energy in joules
- # Define the work distance (height of the electron beam gun)
- # w_d = 12e-3 # Work distance in meters (272 mm)
-
- # Domain setup (Section 2.1.1)
- # Define the 3D domain for the simulation
- x_min, x_max = -0.5 * x_num * dx, 0.5 * x_num * dx # -10 mm to 10 mm
- # print('x_min, x_max:', x_min, x_max)
- y_min, y_max = -0.5 * y_num * dy, 0.5 * y_num * dy # -10 mm to 10 mm
- # print('y_min, y_max:', y_min, y_max)
- z_min, z_max = 0, z_num * dz
- # print('z_min, z_max:', z_min, z_max)
- # x_num = int((x_max - x_min) / dx) + 1
- # y_num = int((y_max - y_min) / dy) + 1
- # z_num = int((z_max - z_min) / dz) + 1
- # print(f"Domain dimensions: x={x_num}, y={y_num}, z={z_num}")
- # domain = np.zeros((z_num, y_num, x_num)) # (z, y , x)
- # print('domain shape: ', np.shape(domain))
- # domain[0:2, :, :] = 1 # Solid cells
- # domain[2, :, :] = 0.5 # Interface layer
- fill_levels_reshaped = np.reshape(fill_levels, (z_num, y_num, x_num))
- # fill_levels_reshaped = fill_levels
- # print('fill_levels_reshaped:', fill_levels_reshaped)
- out_of_range_fill_levels = np.logical_or(fill_levels_reshaped < -tolerance, fill_levels_reshaped > 1 + tolerance)
- if np.any(out_of_range_fill_levels):
- out_of_range_fill_level_values = fill_levels_reshaped[out_of_range_fill_levels]
- out_of_range_fill_level_indices = np.where(out_of_range_fill_levels)
-
- for idx, value in zip(zip(*out_of_range_fill_level_indices), out_of_range_fill_level_values):
- print(f"Index: {idx}, Value: {value}")
- raise ValueError(f"Fill level values must be between 0 and 1 inclusive (with tolerance of {tolerance}).")
-
- sigma = 0.0
- # Different sigma for each dimension
- # sigma = [0.0, 1.0, 0.5] # [z, y, x]
- fill_levels_smoothed = ndimage.gaussian_filter(fill_levels_reshaped, sigma=sigma, mode='nearest', radius=1)
-
- fill_levels_smoothed = np.clip(
- fill_levels_smoothed,
- a_min=np.float64(0.0),
- a_max=np.float64(1.0)
- )
-
- # Check if values are still in bounds
- if not (np.all(fill_levels_smoothed >= 0) and np.all(fill_levels_smoothed <= 1)):
- raise ValueError("Smoothed fill levels out of bounds!")
-
-
- def interpolate_normals(_normal_below: np.ndarray, _normal_above: np.ndarray, t: float) -> np.ndarray:
- """
- Interpolate between two normals using spherical linear interpolation (SLERP)
- """
- # Input validation with better error messages
- if not (0 <= t <= 1):
- raise ValueError(f"Interpolation parameter t must be in [0,1], got {t}")
- if not (isinstance(_normal_below, np.ndarray) and isinstance(_normal_above, np.ndarray)):
- raise ValueError("Normals must be numpy arrays")
-
- # Ensure normals are normalized and correct type
- _normal_below = normalize_vector(_normal_below.astype(np.float64))
- _normal_above = normalize_vector(_normal_above.astype(np.float64))
-
- # Compute dot product with numerical stability
- dot_product = np.clip(np.dot(_normal_below, _normal_above), -1.0, 1.0)
-
- # Better thresholds
- PARALLEL_THRESHOLD = 0.9999 # cos(1°) ≈ 0.9998
- ANTIPARALLEL_THRESHOLD = -0.9999
-
- # Case 1: Normals nearly parallel
- if dot_product > PARALLEL_THRESHOLD:
- _normal_interp = (1.0 - t) * _normal_below + t * _normal_above
- return normalize_vector(_normal_interp)
-
- # Case 2: Normals nearly antiparallel
- if dot_product < ANTIPARALLEL_THRESHOLD:
- # Find perpendicular vector more robustly
- perp = np.array([1.0, 0.0, 0.0])
- if abs(np.dot(_normal_below, perp)) > 0.9:
- perp = np.array([0.0, 1.0, 0.0])
- if abs(np.dot(_normal_below, perp)) > 0.9:
- perp = np.array([0.0, 0.0, 1.0])
- perp = normalize_vector(np.cross(_normal_below, perp))
-
- # Rotate 180 degrees around perp
- angle = np.pi * t
- _normal_interp = _normal_below * np.cos(angle) + perp * np.sin(angle)
- return normalize_vector(_normal_interp)
-
- # Case 3: Standard SLERP with better numerical stability
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
-
- if abs(sin_theta) < 1e-10:
- return normalize_vector(_normal_below)
-
- scale0 = np.sin((1.0 - t) * theta) / sin_theta
- scale1 = np.sin(t * theta) / sin_theta
- _normal_interp = scale0 * _normal_below + scale1 * _normal_above
-
- return normalize_vector(_normal_interp)
-
-
- def interpolate_fill_level_and_normal(_fill_level, surface_normals, _y, _x):
- # Validate input
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
- # Find the indices where the fill level crosses 0.5 in either direction
- indices_above = np.where((_fill_level[:, _y, _x][:-1] > 0.5) & (_fill_level[:, _y, _x][1:] <= 0.5))[0]
- print('indices_above:', indices_above)
- indices_below = np.where((_fill_level[:, _y, _x][:-1] <= 0.5) & (_fill_level[:, _y, _x][1:] > 0.5))[0]
- print('indices_below:', indices_below)
-
- # if len(indices_above) == 0 and len(indices_below) == 0:
- # Fill level does not cross 0.5 at the given (_x, _y) coordinate
- # return None, None, None, None, None, None, None
-
- if len(indices_above) > 0:
- print("len(indices_above) > 0")
- # Fill level crosses 0.5 from above to below
- z_index = indices_above[0]
- print('z_index:', z_index)
- _fill_level_above = _fill_level[z_index, _y, _x]
- print('_fill_level_above:', _fill_level_above)
- _fill_level_below = _fill_level[z_index + 1, _y, _x]
- print('_fill_level_below:', _fill_level_below)
- elif len(indices_below) > 0:
- print("len(indices_below) > 0")
- # Fill level crosses 0.5 from below to above
- z_index = indices_below[0]
- print('z_index:', z_index)
- _fill_level_below = _fill_level[z_index, _y, _x]
- print('_fill_level_below:', _fill_level_below)
- _fill_level_above = _fill_level[z_index + 1, _y, _x]
- print('_fill_level_above:', _fill_level_above)
- else:
- # Fill level does not cross 0.5 at the given (_x, _y) coordinate
- print("No interface cells found")
- return None, None, None, None, None, None, None
-
- # Interpolate the z-coordinate where the fill level is 0.5
- denominator = _fill_level_below - _fill_level_above
- if abs(denominator) < 1e-6:
- # Handle the case where the denominator is very close to zero
- _z_interp = z_index + 0.5
- else:
- _z_interp = z_index + ((0.5 - _fill_level_above) / denominator)
- # _z_interp = z_index + (0.5 - bilinear_interpolate(_x, _y, z_index, _fill_level)) / (bilinear_interpolate(_x, _y, z_index+1, _fill_level) - bilinear_interpolate(_x, _y, z_index, _fill_level))
- print('_z_interp:', _z_interp)
-
- _z_height = (_z_interp + 0.5) * dz
- print('_z_height:', _z_height)
-
- # Interpolate the surface normal at the interpolated z-coordinate
- _normal_below = surface_normals[z_index, _y, _x].astype(np.float64)
- print('_normal_below:', _normal_below)
- _normal_above = surface_normals[z_index + 1, _y, _x].astype(np.float64)
- print('_normal_above:', _normal_above)
-
- t = _z_interp - z_index
-
- # Use spherical interpolation (slerp) for more accurate normal interpolation
- dot_product = np.dot(_normal_below, _normal_above).astype(np.float64)
- if dot_product > np.cos(np.deg2rad(1)): # Threshold for using linear interpolation
- # If normals are very close, use linear interpolation
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- # Use spherical interpolation
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + (np.sin(t*theta) / sin_theta) * _normal_above
- # _normal_interp = quaternion_interpolate(_normal_below, _normal_above, t)
-
- _normal_interp = normalize_vector(_normal_interp)
- print('_normal_interp:', _normal_interp)
-
- return _z_interp, _z_height, _fill_level_below, _fill_level_above, _normal_interp, _normal_below, _normal_above
-
-
- def interpolate_fill_level_and_normal1(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- # Validate input
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at the given (x,y) position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom (reverse order)
- for z in range(len(column)-2, -1, -1): # Start from second-to-last index
- fill_current = column[z]
- fill_next = column[z + 1]
-
- # Check if we cross 0.5 threshold
- if (fill_current - 0.5) * (fill_next - 0.5) <= 0:
- # Found interface
- z_index = z
- _fill_level_above = fill_current
- _fill_level_below = fill_next
-
- print(f"Found interface at z={z_index}")
- print(f"fill_level_above: {_fill_level_above}")
- print(f"fill_level_below: {_fill_level_below}")
-
- # Interpolate z position
- denominator = _fill_level_below - _fill_level_above
- if abs(denominator) < 1e-6:
- _z_interp = z_index + 0.5
- else:
- _z_interp = z_index + ((0.5 - _fill_level_above) / denominator)
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and interpolate normals
- _normal_below = surface_normals[z_index, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z_index + 1, _y, _x].astype(np.float64)
-
- # Interpolate normals
- t = _z_interp - z_index
- dot_product = np.dot(_normal_below, _normal_above)
-
- if dot_product > np.cos(np.deg2rad(1)):
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + \
- (np.sin(t*theta) / sin_theta) * _normal_above
-
- _normal_interp = normalize_vector(_normal_interp)
-
- return _z_interp, _z_height, _fill_level_below, _fill_level_above, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
-
- def interpolate_fill_level_and_normal2(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at the given (x,y) position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom
- for z in range(len(column)-2, -1, -1):
- fill_current = column[z]
- fill_next = column[z + 1]
-
- # Look for significant transition from liquid (>0.5) to gas (<0.2)
- if fill_current > 0.5 and fill_next < 0.2:
- print(f"Found interface at z={z}")
- print(f"fill_level_above: {fill_next}") # gas phase
- print(f"fill_level_below: {fill_current}") # liquid phase
-
- # Interpolate z position
- denominator = fill_next - fill_current
- if abs(denominator) < 1e-6:
- _z_interp = z + 0.5
- else:
- _z_interp = z + ((0.5 - fill_current) / denominator)
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and interpolate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- # Interpolate normal
- t = _z_interp - z
- dot_product = np.dot(_normal_below, _normal_above)
-
- if dot_product > np.cos(np.deg2rad(1)):
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + \
- (np.sin(t*theta) / sin_theta) * _normal_above
-
- _normal_interp = normalize_vector(_normal_interp)
-
- return _z_interp, _z_height, fill_next, fill_current, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
-
- def interpolate_fill_level_and_normal3(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom
- for z in range(len(column)-2, -1, -1):
- current_fill = column[z]
- next_fill = column[z + 1]
-
- # Look for liquid-gas interface (current > 0, next = 0)
- if current_fill > 0 and next_fill == 0:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}") # liquid phase
- print(f"fill_level_above: {next_fill}") # gas phase
-
- # Interface is at z since this is the last liquid cell
- denominator = next_fill - current_fill
- if abs(denominator) < 1e-6:
- _z_interp = z + 0.5
- else:
- _z_interp = z + ((0.5 - current_fill) / denominator)
-
- print('_z_interp', _z_interp)
- _z_height = (_z_interp + 0.5) * dz
- print("_z_height:", _z_height)
-
- # Get and interpolate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- # Interpolate normal
- t = _z_interp - z
-
- # Use default normal if either normal is invalid
- if np.any(np.isnan(_normal_below)) or np.any(np.isnan(_normal_above)):
- _normal_interp = np.array([0.0, 0.0, 1.0])
- else:
- # Interpolate normal at interface
- _normal_interp = interpolate_normals(_normal_below, _normal_above, t)
-
- # _normal_interp = interpolate_normals(_normal_below, _normal_above, t)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
-
- def interpolate_fill_level_and_normal4(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals: {surface_normals.ndim}")
-
- # Get fill levels at position
- column = _fill_level[:, _y, _x]
-
- # Constants
- INTERFACE_THRESHOLD = 0.5 # Standard interface threshold
- GRADIENT_THRESHOLD = 0.1 # Minimum gradient for interface detection
- TOLERANCE = 1e-6 # Numerical tolerance
-
- # Scan from top to bottom
- for z in range(len(column)-2, -1, -1):
- current_fill = column[z]
- next_fill = column[z + 1]
-
- # Check for significant gradient and crossing of 0.5
- gradient = abs(next_fill - current_fill)
- if gradient > GRADIENT_THRESHOLD:
- # Check if we cross the interface threshold
- if ((current_fill > INTERFACE_THRESHOLD > next_fill) or
- (current_fill < INTERFACE_THRESHOLD < next_fill)):
-
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}")
- print(f"fill_level_above: {next_fill}")
-
- # Calculate interpolation parameter
- denominator = next_fill - current_fill
- if abs(denominator) < TOLERANCE:
- t_interface = 0.5
- _z_interp = z + 0.5
- else:
- t_interface = (INTERFACE_THRESHOLD - current_fill) / denominator
- t_interface = np.clip(t_interface, 0.0, 1.0) # Ensure t is in [0,1]
- _z_interp = z + t_interface
-
- print('_z_interp', _z_interp)
- _z_height = (_z_interp + 0.5) * dz
-
- # Get normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- # Handle invalid normals
- if np.any(np.isnan(_normal_below)) or np.any(np.isnan(_normal_above)):
- _normal_interp = np.array([0.0, 0.0, 1.0])
- else:
- # Use clipped t_interface for normal interpolation
- _normal_interp = interpolate_normals(_normal_below, _normal_above, t_interface)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
- def interpolate_fill_level_and_normal5(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals: {surface_normals.ndim}")
-
- # Get fill levels at position
- _column = _fill_level[:, _y, _x]
-
- # Constants with physical meaning
- INTERFACE_THRESHOLD = 0.5 # Standard VOF interface threshold
- GRADIENT_THRESHOLD = 0.1 # Minimum gradient for interface detection
- TOLERANCE = 1e-6 # Numerical tolerance
-
- # Scan from top to bottom
- for z in range(len(_column)-2, -1, -1):
- current_fill = _column[z]
- next_fill = _column[z + 1]
-
- # Check for significant gradient
- gradient = abs(next_fill - current_fill)
- if gradient > GRADIENT_THRESHOLD:
- # Check for interface crossing (VOF method)
- crosses_interface = (
- (current_fill >= INTERFACE_THRESHOLD >= next_fill) or
- (current_fill <= INTERFACE_THRESHOLD <= next_fill)
- )
-
- if crosses_interface:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}")
- print(f"fill_level_above: {next_fill}")
- # Linear interpolation with better numerical stability
- denominator = next_fill - current_fill
- if abs(denominator) < TOLERANCE:
- _z_interp = z + 0.5
- t_interface = 0.5
- else:
- # Calculate interpolation parameter
- t_interface = (INTERFACE_THRESHOLD - current_fill) / denominator
- if not (0 <= t_interface <= 1):
- error_msg = (
- f"\nInterpolation parameter t_interface={t_interface:.6f} is out of bounds [0,1]\n"
- f"Debug information:\n"
- f" - z-index: {z}\n"
- f" - current_fill: {current_fill:.6f}\n"
- f" - next_fill: {next_fill:.6f}\n"
- f" - denominator: {denominator:.6f}\n"
- f" - INTERFACE_THRESHOLD: {INTERFACE_THRESHOLD}\n"
- f" - Position (y,x): ({_y}, {_x})\n"
- f" - Fill level values around interface:\n"
- )
- # Nearby fill levels for context
- for i in range(max(0, z-2), min(len(_column), z+3)):
- error_msg += f" {i:3d} | {_column[i]:.6f}\n"
-
- raise ValueError(error_msg)
- # t_interface = np.clip(t_interface, 0.0, 1.0)
- _z_interp = z + t_interface
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and validate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- _normal_interp = interpolate_normals(_normal_below, _normal_above, t_interface)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- return None, None, None, None, None, None, None
-
-
- def compute_surface_normals_sobel(_fill_levels: np.ndarray) -> np.ndarray:
- """
- Compute surface normals using Sobel filter with enhanced numerical precision.
- """
- # Convert input to high precision at the start
- _fill_levels = _fill_levels.astype(np.float64)
-
- # Input validation
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- # Define tolerances for different checks
- TOLERANCE = {
- 'fill_level_bounds': 1e-10,
- 'zero': 1e-10,
- 'norm': 1e-10,
- 'normalization': 1e-8
- }
-
- # Validate domain bounds with high precision
- '''invalid_elements = np.logical_or(_fill_levels < -TOLERANCE['fill_level_bounds'],
- _fill_levels > 1 + TOLERANCE['fill_level_bounds'])
- if np.any(invalid_elements):
- invalid_indices = np.where(invalid_elements)
- for _z, _y, _x in zip(*invalid_indices):
- _value = _fill_levels[_z, _y, _x]
- print(f"Warning: Invalid value {_value} at [z={_z}, y={_y}, x={_x}]")'''
-
- # Compute gradients using Sobel filter with explicit high precision
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64)
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64)
- _grad_z = -1 * ndimage.sobel(_fill_levels, axis=0, output=np.float64) # multiply by -1 to invert the surface normals to point in +z direction
-
- # Stack gradients into normal vector array with maintained precision
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Compute norm with high precision
- norm = np.linalg.norm(_n_S, axis=-1)
-
- # Create masks for valid regions
- interface_mask = (_fill_levels > TOLERANCE['zero']) & (_fill_levels < 1 - TOLERANCE['zero'])
- nonzero_mask = norm > TOLERANCE['norm']
- valid_mask = interface_mask & nonzero_mask
-
- # Normalize vectors only where valid
- # _n_S[valid_mask] /= norm[valid_mask][..., np.newaxis]
- # _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized = np.full_like(_n_S, np.nan)
- _n_S_normalized[valid_mask] = _n_S[valid_mask] / norm[valid_mask, np.newaxis]
-
- # Verify normalization
- if not np.allclose(np.linalg.norm(_n_S_normalized[valid_mask], axis=-1),1.0, atol=TOLERANCE['normalization']):
- raise ValueError("Normalization failed for some vectors")
-
- # return _n_S
- return _n_S_normalized
-
-
- def compute_surface_normals_sobel1(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'fill_level_bounds': 1e-10,
- 'zero': 1e-10,
- 'norm': 1e-10,
- 'normalization': 1e-8
- }
-
- # Compute gradients
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64)
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64)
- _grad_z = -1 * ndimage.sobel(_fill_levels, axis=0, output=np.float64)
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Compute norm
- norm = np.linalg.norm(_n_S, axis=-1)
-
- # Create masks for valid regions:
- # 1. Interface cells (fill level between 0 and 1)
- interface_mask = (_fill_levels > TOLERANCE['zero']) & (_fill_levels < 1 - TOLERANCE['zero'])
-
- # 2. Cells adjacent to fill level transitions
- # filled_mask = _fill_levels >= 1 - TOLERANCE['zero']
- filled_mask = _fill_levels >= 1
- # empty_mask = _fill_levels <= TOLERANCE['zero']
- empty_mask = _fill_levels <= 0
-
- # Check for transitions in z direction
- z_transitions = np.zeros_like(_fill_levels, dtype=bool)
- z_transitions[:-1, :, :] |= (filled_mask[:-1, :, :] & empty_mask[1:, :, :])
- z_transitions[1:, :, :] |= (empty_mask[:-1, :, :] & filled_mask[1:, :, :])
-
- # Combine masks
- valid_mask = interface_mask | z_transitions
-
- # Default normal for sharp transitions (pointing up)
- _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized[..., 2] = 1.0 # Default z-direction normal
-
- # Expand valid_mask to match _n_S shape
- valid_mask_expanded = valid_mask[..., np.newaxis]
-
- # Where we have valid gradients, use them
- gradient_mask = (norm > TOLERANCE['norm'])[..., np.newaxis]
- mask = valid_mask_expanded & gradient_mask
-
- # Normalize vectors where mask is True
- norm_expanded = norm[..., np.newaxis]
- _n_S_normalized = np.where(mask, _n_S / np.where(norm_expanded > 0, norm_expanded, 1.0), _n_S_normalized)
-
- return _n_S_normalized
-
-
- def compute_surface_normals_sobel2(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'fill_level_bounds': 1e-10,
- 'zero': 1e-10,
- 'norm': 1e-10,
- 'normalization': 1e-8
- }
-
- # Compute gradients
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64)
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64)
- _grad_z = -1 * ndimage.sobel(_fill_levels, axis=0, output=np.float64)
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Compute norm
- norm = np.linalg.norm(_n_S, axis=-1)
-
- # Create masks for valid regions:
- # 1. Interface cells (fill level between 0 and 1)
- interface_mask = (_fill_levels > TOLERANCE['zero']) & (_fill_levels < 1 - TOLERANCE['zero'])
-
- # 2. Cells adjacent to fill level transitions
- filled_mask = _fill_levels >= 1 - TOLERANCE['zero']
- empty_mask = _fill_levels <= TOLERANCE['zero']
-
- # Check for transitions in z direction
- z_transitions = np.zeros_like(_fill_levels, dtype=bool)
- z_transitions[:-1, :, :] |= (filled_mask[:-1, :, :] & empty_mask[1:, :, :])
- z_transitions[1:, :, :] |= (empty_mask[:-1, :, :] & filled_mask[1:, :, :])
-
- # Combine masks
- valid_mask = interface_mask | z_transitions
-
- # Initialize output array with NaN
- _n_S_normalized = np.full_like(_n_S, np.nan)
-
- # Create mask for valid gradients
- gradient_mask = norm > TOLERANCE['norm']
-
- # Combine masks and reshape for broadcasting
- mask = valid_mask & gradient_mask
-
- # Normalize only where mask is True
- for i in range(3): # For each component (x, y, z)
- _n_S_normalized[..., i] = np.where(mask,
- _n_S[..., i] / np.where(norm > 0, norm, 1.0),
- np.nan)
-
- return _n_S_normalized
-
- def compute_surface_normals_sobel3(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with default upward normal
- _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized[..., 2] = 1.0
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
-
- def interpolate_normals_new(_normal_below: np.ndarray, _normal_above: np.ndarray, t: float) -> np.ndarray:
- """Interpolate between two normals using SLERP"""
- if not (0 <= t <= 1):
- raise ValueError(f"Interpolation parameter t must be in [0,1], got {t}")
-
- # Ensure normalized vectors
- _normal_below = normalize_vector(_normal_below.astype(np.float64))
- _normal_above = normalize_vector(_normal_above.astype(np.float64))
-
- # Compute dot product
- dot_product = np.clip(np.dot(_normal_below, _normal_above), -1.0, 1.0)
-
- # Thresholds for special cases
- PARALLEL_THRESHOLD = 0.9999
- ANTIPARALLEL_THRESHOLD = -0.9999
-
- # Case 1: Nearly parallel
- if dot_product > PARALLEL_THRESHOLD:
- return normalize_vector((1.0 - t) * _normal_below + t * _normal_above)
-
- # Case 2: Nearly antiparallel
- if dot_product < ANTIPARALLEL_THRESHOLD:
- perp = np.array([1.0, 0.0, 0.0])
- if abs(np.dot(_normal_below, perp)) > 0.9:
- perp = np.array([0.0, 1.0, 0.0])
- perp = normalize_vector(np.cross(_normal_below, perp))
- angle = np.pi * t
- return normalize_vector(_normal_below * np.cos(angle) + perp * np.sin(angle))
-
- # Case 3: Standard SLERP
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
-
- if abs(sin_theta) < 1e-10:
- return normalize_vector(_normal_below)
-
- scale0 = np.sin((1.0 - t) * theta) / sin_theta
- scale1 = np.sin(t * theta) / sin_theta
- return normalize_vector(scale0 * _normal_below + scale1 * _normal_above)
-
-
- def interpolate_fill_level_and_normal_new(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals: {surface_normals.ndim}")
-
- # Get fill levels at position
- _column = _fill_level[:, _y, _x]
-
- # Constants with physical meaning
- INTERFACE_THRESHOLD = 0.5 # Standard VOF interface threshold
- GRADIENT_THRESHOLD = 0.1 # Minimum gradient for interface detection
- TOLERANCE = 1e-6 # Numerical tolerance
-
- # Scan from top to bottom
- for z in range(len(_column)-2, -1, -1):
- current_fill = _column[z]
- next_fill = _column[z + 1]
-
- # Check for significant gradient
- gradient = abs(next_fill - current_fill)
- if gradient > GRADIENT_THRESHOLD:
- # Check for interface crossing (VOF method)
- '''crosses_interface = (
- (current_fill >= INTERFACE_THRESHOLD - TOLERANCE >= next_fill) or
- (current_fill <= INTERFACE_THRESHOLD + TOLERANCE <= next_fill)
- )'''
- crosses_interface = (
- (current_fill >= (INTERFACE_THRESHOLD - TOLERANCE) and
- next_fill <= (INTERFACE_THRESHOLD + TOLERANCE))
- )
-
- if crosses_interface:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}")
- print(f"fill_level_above: {next_fill}")
- # Linear interpolation with better numerical stability
- denominator = next_fill - current_fill
- if abs(denominator) < TOLERANCE:
- _z_interp = z + 0.5
- t_interface = 0.5
- else:
- # Calculate interpolation parameter
- t_interface = (INTERFACE_THRESHOLD - current_fill) / denominator
- if not (0 <= t_interface <= 1):
- error_msg = (
- f"\nInterpolation parameter t_interface={t_interface:.6f} is out of bounds [0,1]\n"
- f"Debug information:\n"
- f" - z-index: {z}\n"
- f" - current_fill: {current_fill:.6f}\n"
- f" - next_fill: {next_fill:.6f}\n"
- f" - denominator: {denominator:.6f}\n"
- f" - INTERFACE_THRESHOLD: {INTERFACE_THRESHOLD}\n"
- f" - Position (y,x): ({_y}, {_x})\n"
- f" - Fill level values around interface:\n"
- )
- # Nearby fill levels for context
- for i in range(max(0, z-2), min(len(_column), z+3)):
- error_msg += f" {i:3d} | {_column[i]:.6f}\n"
-
- raise ValueError(error_msg)
- # t_interface = np.clip(t_interface, 0.0, 1.0)
- _z_interp = z + t_interface
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and validate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- _normal_interp = interpolate_normals_new(_normal_below, _normal_above, t_interface)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- return None, None, None, None, None, None, None
-
-
- def compute_surface_normals_sobel_zeros(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with default upward normal
- _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized[..., 2] = 1.0
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
-
- def compute_surface_normals_sobel_nans(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter with NaN initialization"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with NaN
- _n_S_normalized = np.full_like(_n_S, np.nan)
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
- # Compute surface normals (Section 2.1.1)
- # n_S, (grad_x, grad_y, grad_z) = compute_surface_normals(domain)
- n_S = compute_surface_normals_sobel_nans(fill_levels_smoothed)
- # print('n_S: ', np.shape(n_S))
-
- # Beam profile (Section 2.1)
- # Define the Gaussian beam profile
- sigma: float = 500e-6 / (4 * dx) # Beam profile sigma
- x = np.arange(-np.floor(min(x_num//4, 2*sigma)), 1+np.floor(min((x_num//4), 2*sigma)))
- # print('x:', x)
- y = np.arange(-floor(min(y_num//4, 2*sigma)), 1+floor(min((y_num//4), 2*sigma)))
- # print('y:', y)
- X, Y = np.meshgrid(x, y)
- beam_profile = np.exp(-(X ** 2 + Y ** 2) / (2 * sigma ** 2))
- beam_profile /= np.sum(beam_profile) # Normalize beam profile
- print('beam_profile shape: ', np.shape(beam_profile))
-
- beam_size = np.shape(beam_profile)[0] * np.shape(beam_profile)[1]
- # print('Beam size:', beam_size)
-
- # Detector setup (Section 2.1.1)
- # Define positions of the detectors
- detector_height = 4 * z_max # Detector height in meters (272 mm)
- detector_spacing_x = 1.5 * dx * x_num # Detector spacing in meters (127 mm)
- detector_spacing_y = detector_spacing_x # 1.5 * dy * y_num # Detector spacing in meters (127 mm)
- detector_diameter = 4 * dx # Detector diameter in meters (50 mm)
-
- # Define detector positions relative to the origin
- detector_positions = np.array([
- [0.5 * detector_spacing_x, 0, detector_height],
- [-0.5 * detector_spacing_x, 0, detector_height],
- [0, 0.5 * detector_spacing_y, detector_height],
- [0, -0.5 * detector_spacing_y, detector_height]
- ], dtype=np.float64)
- # print('detector_positions:', detector_positions)
-
-
- # Define detector vertices for solid angle calculation
- def get_detector_vertices(detector_position: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
- # Calculate the side length of the detector
- side_length : float = np.sqrt(3) * (detector_diameter / 2)
- height : float = np.sqrt(3) * (side_length / 2)
-
- # Determine the orientation based on the detector position
- if np.allclose(detector_position[:2], [0, 0]):
- raise ValueError("Detector position cannot be at the origin (0, 0, 0)")
- elif detector_position[0] > 0: # Right detector
- _v1 = detector_position + np.array([2 / 3 * height, 0, 0], dtype=np.float64)
- _v2 = detector_position + np.array([-height / 3, side_length / 2, 0], dtype=np.float64)
- _v3 = detector_position + np.array([-height / 3, -side_length / 2, 0], dtype=np.float64)
- elif detector_position[0] < 0: # Left detector
- _v1 = detector_position + np.array([-2 / 3 * height, 0, 0], dtype=np.float64)
- _v2 = detector_position + np.array([height / 3, -side_length / 2, 0], dtype=np.float64)
- _v3 = detector_position + np.array([height / 3, side_length / 2, 0], dtype=np.float64)
- elif detector_position[1] > 0: # Back detector
- _v1 = detector_position + np.array([0, 2 / 3 * height, 0], dtype=np.float64)
- _v2 = detector_position + np.array([-side_length / 2, -height / 3, 0], dtype=np.float64)
- _v3 = detector_position + np.array([side_length / 2, -height / 3, 0], dtype=np.float64)
- elif detector_position[1] < 0: # Front detector
- _v1 = detector_position + np.array([0, -2 / 3 * height, 0], dtype=np.float64)
- _v2 = detector_position + np.array([side_length / 2, height / 3, 0], dtype=np.float64)
- _v3 = detector_position + np.array([-side_length / 2, height / 3, 0], dtype=np.float64)
- else:
- raise ValueError("Invalid detector position. Must be in one of the four expected positions.")
-
- return _v1, _v2, _v3
-
-
- # Compute solid angle using Equation A1 (Appendix A)
- def compute_solid_angle(vector1: np.ndarray, vector2: np.ndarray, vector3: np.ndarray) -> float:
- # Input validation
- if not isinstance(vector1, np.ndarray) or not isinstance(vector2, np.ndarray) or not isinstance(vector3, np.ndarray):
- raise ValueError("Inputs must be numpy arrays")
- if vector1.shape != (3,) or vector2.shape != (3,) or vector3.shape != (3,):
- raise ValueError("Inputs must be 3D vectors")
- if not np.isfinite(vector1).all() or not np.isfinite(vector2).all() or not np.isfinite(vector3).all():
- raise ValueError("Inputs must be finite")
-
- # Convert to high precision once
- vectors = [vec.astype(np.float64) for vec in (vector1, vector2, vector3)]
- cross_product = np.cross(vectors[1], vectors[2])
- # Check if the cross product is nearly zero
- if np.linalg.norm(cross_product) < 1e-16:
- raise ValueError("Vectors are nearly coplanar")
-
- # Compute solid angle
- try:
- # Compute the numerator of the solid angle formula
- numerator: float = np.dot(vectors[0], cross_product)
- # Check if the numerator is nearly zero
- if np.abs(numerator) < 1e-16:
- raise ValueError("Vectors are nearly parallel")
-
- norms = [np.linalg.norm(v) for v in vectors]
-
- # Compute the denominator of the solid angle formula
- denominator = (
- norms[0] * norms[1] * norms[2] +
- np.dot(vectors[0], vectors[1]) * norms[2] +
- np.dot(vectors[0], vectors[2]) * norms[1] +
- np.dot(vectors[1], vectors[2]) * norms[0]
- )
-
- # Compute the solid angle in radians
- solid_angle = 2 * np.arctan2(numerator, denominator)
- # print("Solid angle:", solid_angle)
- return solid_angle
- except ZeroDivisionError:
- raise ValueError("Error computing solid angle: division by zero")
- except np.linalg.LinAlgError:
- raise ValueError("Error computing solid angle: linear algebra error")
-
-
- # BSE coefficient functions (Equations 11 and 10, Section 2.1.2)
- def eta_0(_A: float) -> float:
- """Equation 11 from the paper"""
- _eta_0: float = -0.0254 + 0.016 * _A - 1.86e-4 * _A ** 2 + 8.3e-7 * _A ** 3
- print('eta_0: ', _eta_0)
- return _eta_0
-
-
- # BSE Coefficient for single atomic targets, found to hold in the range of 10–100 keV
- def eta(theta: float, _A: float) -> float:
- """Equation 10 from the paper"""
- B: float = 0.89
- _eta: float = B * (eta_0(_A) / B) ** np.cos(theta)
- print('eta: ', _eta)
- return _eta
-
-
- '''
- # Compute the weighted average BSE coefficient for the alloy
- def compute_weighted_eta(theta: float, alloy: Alloy) -> float:
- """Equation 12 from the paper"""
- weighted_eta: float = 0
- for element, weight_fraction in zip(alloy.elements, alloy.composition):
- # print('element: ', element)
- # print('weight_fraction: ', weight_fraction)
- _A: float = element.atomic_number
- # print('A: ', _A)
- eta_i: float = eta(theta, _A)
- # print('eta_i: ', eta_i)
- weighted_eta += weight_fraction * eta_i
- # print('weighted_eta: ', weighted_eta)
- return weighted_eta
- '''
-
-
- # Correction function C(theta) (Equation 14, Section 2.1.3)
- def C(theta: float) -> float:
- print("theta1:", theta)
- """Equation 14 from the paper"""
- theta = np.rad2deg(theta) # Ensure angle is in degrees for C function
- _C: float = 1.01 - 4.06 * theta + 1.36e-5 * theta ** 2 - 1.71e-6 * theta ** 3
- print('C: ', _C)
- # Check if C is negative and raise an exception
- if _C < 0:
- raise ValueError(f"C function returned a negative value for theta: {theta}. "
- f"Please check the input and the C function implementation.")
- return _C
-
-
- # Scattering function g_BSE (Equation 13, Section 2.1.3)
- def g_BSE(_theta: float, _alpha: float, _beta: float) -> float:
- print("_theta, _alpha, _beta:", _theta, _alpha, _beta)
- """Equation 13 from the paper"""
- # k: float = np.rad2deg(_theta) / 13 # Correct implementation as per the paper
- # print("k: ", k)
- diffusive_part: float = (eta_0(A) / pi) * np.cos(_alpha)
- print(np.cos(_alpha))
- print("diffusive_part:", diffusive_part)
- # reflective_part: float = ((eta(_theta, A) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
- # if np.abs(reflective_part) < 1e-10:
- # reflective_part = 0
- # print("reflective_part:", reflective_part)
-
- # Compute a weighted average BSE coefficient for the alloy
- # reflective_part: float = ((compute_weighted_eta(_theta, Ti6Al4V) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
-
- # print('g_BSE: ', diffusive_part + reflective_part)
- return diffusive_part #+ reflective_part
-
-
- # Calculate the beam direction vector p_E
- def calculate_p_E(_beam_x, _beam_y):
- print(f"calculate_p_E({_beam_x}, {_beam_y}):")
- # Calculate the vector from the beam gun to the beam location
- beam_vector = np.array([_beam_x * dx, _beam_y * dy, -detector_height], dtype=np.float64)
- return normalize_vector(beam_vector)
-
- def calculate_p_E1(_beam_x, _beam_y, _height):
- """Calculate normalized beam direction vector from gun to surface point
-
- Args:
- _beam_x: Beam x-position relative to center (in cells)
- _beam_y: Beam y-position relative to center (in cells)
- _height: Height difference between gun and surface point (m)
- """
- print(f"calculate_p_E({_beam_x}, {_beam_y}):")
- print("_height:", _height)
-
- # Calculate vector components
- x_component = _beam_x * dx # Convert cells to meters
- y_component = _beam_y * dy # Convert cells to meters
- z_component = _height # Already in meters
-
- # Create beam vector (from gun to surface point)
- beam_vector = np.array([x_component, y_component, z_component], dtype=np.float64)
-
- # Verify downward direction
- if normalize_vector(beam_vector)[2] > 0:
- raise ValueError("Beam vector pointing upward")
-
- # Normalize and return
- return normalize_vector(beam_vector)
-
-
- # Calculate phi as the angle between p_E and the negative z-axis
- def calculate_phi(_p_E) -> float:
- z_axis = np.array([0.0, 0.0, -1.0], dtype=np.float64)
-
- # Check if _p_E is approximately equal to the negative z-axis
- '''if np.allclose(_p_E, z_axis, rtol=1e-8, atol=1e-8):
- _phi_rad = 0.0
- else:
- # Calculate the angle between _p_E and the z-axis using arccos
- _phi_rad = calculate_angle(_p_E, z_axis)'''
- _phi_rad = calculate_angle(_p_E, z_axis)
- print(f"calculate_phi({_p_E}, {_phi_rad}):")
- return _phi_rad
-
-
- def set_beam_movement(_x_center, _y_center, _x_offset, _y_offset, _movement_type):
- if _movement_type == 'X':
- if _x_offset >= 0:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset - 1
- else:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset + 1
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center + _y_offset + 1
- elif _movement_type == 'Y':
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center + _x_offset + 1
- if _y_offset >= 0:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset - 1
- else:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset + 1
- elif _movement_type == 'XY':
- if _x_offset >= 0:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset - 1
- else:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset + 1
- if _y_offset >= 0:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset - 1
- else:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset + 1
- else:
- raise ValueError("movement_type must be 'X' or 'Y' or 'XY'")
-
- return _beam_x_start, _beam_x_end, _beam_y_start, _beam_y_end
-
-
- # Initialize arrays to store electron counts
- electrons_per_detector_per_step = np.zeros((x_num, y_num, len(detector_positions)))
- total_electrons_per_step = np.zeros((x_num, y_num))
-
- # Additional arrays to capture detailed information
- electrons_per_detector_per_beam_profile = np.zeros(
- (beam_profile.shape[0], beam_profile.shape[1], len(detector_positions)))
- total_electrons_per_beam_profile = np.zeros((beam_profile.shape[0], beam_profile.shape[1]))
-
- # Calculate the center indices
- x_center = x_num // 2
- y_center = y_num // 2
- # print('x_center, y_center:', x_center, y_center)
-
- x_offset: int = 0
- y_offset: int = 0
- movement_type: str = 'Y'
-
- # Beam X-Movement
- # beam_x_start, beam_x_end = x_center - x_offset, x_center + x_offset + 1
- # beam_y_start, beam_y_end = y_center + y_offset, y_center + y_offset + 1
-
- # Beam Y-Movement
- # beam_x_start, beam_x_end = x_center + x_offset, x_center + x_offset + 1
- # beam_y_start, beam_y_end = y_center - y_offset, y_center + y_offset + 1
-
- # For X or Y Movement
- beam_x_start, beam_x_end, beam_y_start, beam_y_end = set_beam_movement(
- x_center, y_center, x_offset, y_offset, movement_type)
-
- # Determine the step size based on the direction of movement
- step_x = 1 if beam_x_start <= beam_x_end else -1
- step_y = 1 if beam_y_start <= beam_y_end else -1
-
- # print(beam_x_start, beam_x_end, beam_y_start, beam_y_end, step_x, step_y)
-
-
- z_heights = np.zeros(beam_size)
- fill_level_b = np.zeros(beam_size)
- fill_level_a = np.zeros(beam_size)
- normals = np.zeros((beam_size, 3))
- normals_b = np.zeros((beam_size, 3))
- normals_a = np.zeros((beam_size, 3))
- cell_counter = 0
-
-
- def normalize_vector(vector: np.ndarray) -> np.ndarray:
- """Normalize a 3D vector."""
- vector = vector.astype(np.float64) # Convert to float64
- norm = np.linalg.norm(vector)
- if norm < 1e-18: # Avoid division by zero
- return vector # Return original if nearly zero
- return vector / norm # Normalize
-
-
- def calculate_angle(_v1: np.ndarray, _v2: np.ndarray) -> float:
- """Calculate the angle between two normalized vectors."""
- _v1 = normalize_vector(_v1)
- _v2 = normalize_vector(_v2)
-
- dot_product = np.clip(np.dot(_v1, _v2), -1.0, 1.0)
-
- return np.arccos(dot_product) # Return angle in radians
-
-
- # Main simulation loop to iterate over the simulation domain
- for beam_x in range(beam_x_start, beam_x_end, step_x):
- for beam_y in range(beam_y_start, beam_y_end, step_y):
- beam_loc = (beam_x, beam_y)
- print(f"Beam center is at domain location: {beam_loc}")
-
- # Calculate the beam direction vector p_E
- # >> p_E = calculate_p_E(beam_x - x_center, beam_y - y_center)
- # >> print('p_E: ', p_E)
-
- # Calculate phi as the angle between p_E and the negative z-axis
- # phi_rad = np.arctan2(np.sqrt(dx**2 + dy**2), w_d)
- # >> phi_rad = calculate_phi(p_E)
- # print(f"phi_rad: {phi_rad}")
-
- # Reset counters for this beam location
- total_electrons: int = 0
- electrons_per_detector = np.zeros(len(detector_positions))
-
- # Iterate over the beam profile
- # for profile_x in [161]:
- for profile_x in range(beam_profile.shape[0]):
- # for profile_y in [42]:
- for profile_y in range(beam_profile.shape[1]):
- print(f"Beam profile shape ({profile_x}, {profile_y})")
- domain_x = beam_loc[0] - beam_profile.shape[0] // 2 + profile_x
- domain_y = beam_loc[1] - beam_profile.shape[1] // 2 + profile_y
- print(f"Domain ({domain_x}, {domain_y})")
-
- # p_E = calculate_p_E(domain_x - x_center, domain_y - y_center)
- # print('p_E: ', p_E)
- # phi_rad = calculate_phi(p_E)
-
- if (0 <= domain_x < x_num) and (0 <= domain_y < y_num):
- if beam_profile[profile_x, profile_y] > threshold:
- print('finding interface_z_indices')
- # print('interface_z_indices:', interface_z)
- z_interp, z_height, fill_level_below, fill_level_above, normal_interp, normals_below, normals_above \
- = interpolate_fill_level_and_normal_new(fill_levels_smoothed, n_S, domain_y, domain_x)
- if z_interp is not None:
- z_heights[cell_counter] = z_height
- fill_level_b[cell_counter] = fill_level_below
- fill_level_a[cell_counter] = fill_level_above
- normals[cell_counter] = normal_interp
- normals_b[cell_counter] = normals_below
- normals_a[cell_counter] = normals_above
- cell_counter += 1
- n_S_local = normal_interp
- print('n_S_local:', n_S_local)
- if not np.allclose(np.linalg.norm(n_S_local), 1.0, atol=1e-5):
- raise ValueError(f"n_S_local is not normalized, {np.linalg.norm(n_S_local)}")
-
- p_E = calculate_p_E1(domain_x - x_center, domain_y - y_center, -(detector_height-z_height))
- print('p_E: ', p_E)
-
- theta_local = calculate_angle(-p_E, n_S_local)
- print('theta_local:', theta_local)
- # if theta_local < 0.0 or theta_local > 0.011:
- # raise ValueError(f"theta_local value should be 0.0, found {theta_local}")
- # theta_local: 0.016
-
- scattering_point = np.array([
- (domain_x - beam_x) * dx,
- (domain_y - beam_y) * dy,
- z_height
- ], dtype=np.float64)
- print(f'scattering_point ({domain_x}, {domain_y}, {z_interp}):', scattering_point)
-
- for det_idx, det_pos in enumerate(detector_positions):
- # print(f"Detector {det_idx + 1} position: {det_pos}")
- det_idx: int
- # det_pos: np.ndarray[float]
- print(f'det_pos[{det_idx + 1}]: {det_pos}')
-
- # Calculate the vector from the scattering point to the detector
- det_vec = det_pos - scattering_point
- det_vec = normalize_vector(det_vec)
- print(f'det_vec[{det_idx + 1}]: {det_vec}')
- if not np.allclose(np.linalg.norm(det_vec), 1.0):
- raise ValueError("det_vec is not normalized.")
-
- # alpha = np.arccos(np.clip(np.dot(det_vec, n_S_local), -1.0, 1.0))
- alpha = calculate_angle(det_vec, n_S_local)
- print(f'alpha: {alpha}')
-
- if alpha > np.pi / 2:
- continue
-
- p_E_reflected = p_E - 2.0 * np.dot(p_E, n_S_local) * n_S_local
- # p_E_reflected = np.where(np.abs(p_E_reflected) < 1e-8, 0.0, p_E_reflected)
- p_E_reflected = normalize_vector(p_E_reflected)
- print(f'p_E_reflected: {p_E_reflected}')
-
- # beta = np.arccos(np.clip(np.dot(p_E_reflected, det_vec), -1.0, 1.0))
- beta = calculate_angle(p_E_reflected, det_vec)
- print(f'beta: {beta}')
-
- # Apply energy threshold
- # electron_energy = N_0 * e * g_bse_value
- # print('electron_energy: ', electron_energy)
-
- # Calculate the initial energy of the primary electrons (E_0)
- V_acc = 60e3 # Accelerating voltage in Volts
- E_0 = V_acc * e # Initial energy of primary electrons in joules
-
- # Calculate the energy of the backscatter electrons
- E_BSE = E_0 * eta(theta_local, A)
- # Compute a weighted average BSE coefficient for the alloy
- # E_BSE = E_0 * compute_weighted_eta(theta_local, Ti6Al4V)
- # print('E_BSE: ', E_BSE)
-
- # Bifurcation logic for PE and SE
- # if E_BSE > threshold_energy:
- # Backscatter electrons (BSE)
- # Compute scattering function g_bse_value
- g_bse_value = g_BSE(theta_local, alpha, beta)
- print(f'g_bse_value: {g_bse_value}')
-
- # Ensure g_bse_value is non-negative
- if g_bse_value <= 0:
- raise ValueError(f"g_bse_value is {g_bse_value} (zero negative), check calculations.")
-
- # Compute solid angle (Equation A1)
- v1, v2, v3 = get_detector_vertices(det_pos)
- normal = np.cross(v2 - v1, v3 - v1)
- normal = normalize_vector(normal)
-
- dOmega = compute_solid_angle(v1 - scattering_point,
- v2 - scattering_point,
- v3 - scattering_point)
- print(f"Detector {det_idx + 1} solid angle: {dOmega}")
-
- # Compute number of electrons collected by the detector (Equation A2)
- # The beam ray stays on each grid point ij for the dwell time t_d and distributes N_0 PE into it
- N_0 = beam_current * dwell_time / e # Equation 3, Section 2.1.1
- print("N_0:", N_0)
-
- # There are N1 electrons in the chamber after the 1st scattering event
- N_1 = N_0 * eta(theta_local, A) # Equation 4, Section 2.1.1
- # Compute a weighted average BSE coefficient for the alloy
- # N_1 = N_0 * compute_weighted_eta(theta_local, Ti6Al4V) # Equation 4, Section 2.1.1
- print("N_1:", N_1)
-
- # The number of electrons N_1_d collected by a general detector d
- N_1_d = N_1 * np.sum(g_bse_value * dOmega) # Equation 5, Section 2.1.1
-
- # Accumulate electrons per detector for this beam profile point
- electrons_per_detector_per_beam_profile[profile_x, profile_y, det_idx] += N_1_d
-
- # Accumulate electrons per detector for this domain point
- electrons_per_detector_per_step[beam_x, beam_y, det_idx] += N_1_d
-
- # Add to total electrons for this beam profile point
- total_electrons_per_beam_profile[profile_x, profile_y] += N_1_d
-
- # Add to total electrons for this domain point
- total_electrons_per_step[beam_x, beam_y] += N_1_d
-
- # Accumulate electrons per detector
- electrons_per_detector[det_idx] += N_1_d
-
- # Add to total electrons
- total_electrons += N_1_d
-
- # print(f"After addition: electrons_per_detector[{det_idx + 1}]: {electrons_per_detector[det_idx]}")
- else:
- print(f"No interface found at ({domain_x}, {domain_y})")
-
- # Print detailed information for this beam profile point
- # print(f"Electrons per detector at this beam profile point ({profile_x}, {profile_y}): {electrons_per_detector_per_beam_profile[profile_x, profile_y]}")
- # print(f"Total electrons at this beam profile point ({profile_x}, {profile_y}): {total_electrons_per_beam_profile[profile_x, profile_y]}")
-
- # Print detailed information for this domain point after iterating over the beam profile
- # print(f"Electrons per detector at this domain point ({beam_x}, {beam_y}): {electrons_per_detector_per_step[beam_x, beam_y]}")
- # print(f"Total electrons at this domain point ({beam_x}, {beam_y}): {total_electrons_per_step[beam_x, beam_y]}")
-
- # print("Shape of electrons_per_detector_per_step:", electrons_per_detector_per_step.shape)
- # print("Number of lines:", len(lines))
-
- # Compute final results
- total_electrons_all_steps = np.sum(total_electrons_per_step)
- total_electrons_per_detector = np.sum(electrons_per_detector_per_step, axis=(0, 1))
-
- # Print final results
- # print(f"Total electrons hitting all detectors: {total_electrons_all_steps}")
- # for i, electrons in enumerate(total_electrons_per_detector):
- # print(f"Total electrons hitting detector {i + 1}: {electrons}")
-
- # for beam_x in range(beam_x_start, beam_x_end, step_x):
- # for beam_y in range(beam_y_start, beam_y_end, step_y):
- # print(f"Electrons per detector at this domain point ({beam_x}, {beam_y}): {electrons_per_detector_per_step[beam_x, beam_y]}")
- # print('end simulation')
- # print(electrons_per_detector_per_step)
-
- # print("total_electrons:", total_electrons, "at timestep", current_timestep)
- # print("electrons_per_det:", electrons_per_det, "at timestep", current_timestep)
- # print("electrons_per_beam_loc:", electrons_per_beam_loc, "at timestep", current_timestep)
- # print("electrons_per_beam_loc_per_det:", electrons_per_beam_loc_per_det, "at timestep", current_timestep)
- return electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, normals, normals_b, normals_a
-
-def read_fill_levels(filepath):
- """Read only the FillingFrac grid from VDB file"""
- try:
- # Try to read just the FillingFrac grid
- grid = vdb.read(filepath, "FillingFrac")
-
- # Get dimensions
- dims = grid.evalActiveVoxelDim()
- print(f"Grid dimensions: {dims}")
-
- # Create array
- fill_array = np.zeros(dims, dtype=np.float64)
-
- # Copy data
- grid.copyToArray(fill_array)
-
- # Transpose to ray tracer's expected order (z, y, x)
- fill_array = np.transpose(fill_array, (2, 1, 0))
-
- # Create standardized array with 256 z cells
- standardized_array = np.zeros((256, dims[1], dims[0]), dtype=np.float64)
-
- # Copy data from original array to standardized array
- standardized_array[:fill_array.shape[0], :, :] = fill_array
-
- # Apply Gaussian smoothing
- sigma = 1.0
- standardized_array_smoothed = ndimage.gaussian_filter(standardized_array, sigma=sigma, mode='nearest', radius=1)
-
- # Print values
- x = 157
- y = 196
- print(f"\nFill levels at ({x}, {y}):")
- print("z_index | fill_level")
- print("-" * 20)
-
- # Get fill levels
- fill_levels = standardized_array_smoothed[:, y, x]
- valid_levels = fill_levels[(fill_levels > 0) & (fill_levels < 1)]
-
- if len(valid_levels) > 0:
- for z, val in enumerate(fill_levels):
- if 0 <= val <= 1:
- print(f"{z:7d} | {val:.6f}")
- else:
- print("No fill level found between 0 and 1")
-
- return standardized_array_smoothed
-
- except Exception as e:
- print(f"Error reading {filepath}")
- print(f"Error details: {str(e)}")
- try:
- metadata = vdb.readMetadata(filepath)
- print("\nFile metadata:")
- for key, value in metadata.items():
- print(f"{key}: {value}")
- except Exception as meta_e:
- print(f"Error reading metadata: {str(meta_e)}")
- return None
-
-def process_single_timestep(filepath, timestep):
- """Process a single VDB file and run ray tracer"""
- try:
- # Read the FillingFrac grid
- grid = vdb.read(filepath, "FillingFrac")
-
- # Get dimensions
- dims = grid.evalActiveVoxelDim()
- x_num, y_num = dims[0], dims[1]
- z_num = 256 # dims[2] # Fixed z dimension
- print(f"Domain dimensions: {x_num} x {y_num} x {z_num}")
-
- # Read fill levels
- fill_levels = read_fill_levels(filepath)
- if fill_levels is None:
- raise ValueError("Could not read fill levels from VDB file")
-
- if fill_levels is not None:
- print("\nSuccessfully read fill levels")
- print(f"Array shape: {np.shape(fill_levels)}")
-
- try:
- # Convert to numpy array if not already
- fill_array = np.asarray(fill_levels)
-
- # Get non-zero elements
- non_zero = fill_array[fill_array > 0]
-
- # Get interface cells (0 < fill_level < 1)
- interface_cells = fill_array[(fill_array > 0) & (fill_array < 1)]
-
- # Print statistics
- if len(non_zero) > 0:
- print("\nNon-zero cells statistics:")
- print(f"Number of non-zero cells: {len(non_zero)}")
- print(f"Min non-zero value: {np.min(non_zero):.6f}")
- print(f"Max value: {np.max(fill_array):.6f}")
-
- if len(interface_cells) > 0:
- print("\nInterface cells statistics:")
- interface_percentage = (len(interface_cells) * 100.0) / len(non_zero)
- print(f"Number of interface cells: {len(interface_cells)} ({interface_percentage:.3f}%)")
- print(f"Min interface value: {np.min(interface_cells):.6f}")
- print(f"Max interface value: {np.max(interface_cells):.6f}")
- print(f"Mean interface value: {np.mean(interface_cells):.6f}")
- else:
- print("\nNo interface cells found (no fill levels between 0 and 1)")
-
- except Exception as e:
- print(f"Error processing array: {str(e)}")
-
- # Run ray tracer
- print(f"\nProcessing timestep {timestep}")
- electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, \
- normals, normals_b, normals_a = runRayTracer(
- x_num,
- y_num,
- z_num,
- # fill_levels.flatten(),
- fill_levels,
- timestep
- )
-
- return electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, \
- normals, normals_b, normals_a
-
- except Exception as e:
- print(f"Error processing timestep {timestep}: {str(e)}")
- raise
-
-def main():
- """Main function to process VDB files and run ray tracer"""
- # Directory containing VDB files
- directory = "/local/ca00xebo/softwares/KiSSAM/Markl-plate"
-
- # Get sorted list of VDB files
- vdb_files = []
- for f in os.listdir(directory):
- if f.startswith('geometry-000003') and f.endswith('.vdb'):
- try:
- # Extract timestep number from filename
- timestep = int(f.split('-')[1].split('.')[0])
- vdb_files.append((timestep, f))
- except (IndexError, ValueError) as e:
- print(f"Skipping file {f}: {str(e)}")
-
- # Sort by timestep
- vdb_files.sort() # sort based on timestep number
-
- total_timesteps = len(vdb_files)
- print(f"Found {total_timesteps} VDB files to process")
- electrons_across_timesteps = np.zeros((total_timesteps, 4))
-
- # Process each timestep
- _results = []
- for timestep, vdb_file in enumerate(vdb_files):
- filepath = os.path.join(directory, vdb_file[1]) # vdb_file[1] contains the filename
- try:
- print(f"\nProcessing file {vdb_file[1]} (timestep {timestep})")
-
- # Process this timestep
- _result = process_single_timestep(filepath, timestep)
-
- # Debug electron counts
- print("\nElectron count debug:")
- print(f"Shape of electrons_per_detector_per_step: {_result[0].shape}")
- print(f"Min value: {np.min(_result[0])}")
- print(f"Max value: {np.max(_result[0])}")
- print(f"Mean value: {np.mean(_result[0])}")
-
- # Store electron counts with debugging
- total_electrons_for_detectors = np.sum(_result[0], axis=(0, 1))
- print(f"Total electrons before assignment: {total_electrons_for_detectors}")
-
- electrons_across_timesteps[timestep] = total_electrons_for_detectors
- print(f"Electrons for timestep {timestep}:")
- print(f"Raw values: {electrons_across_timesteps[timestep]}")
- print(f"Sum: {np.sum(electrons_across_timesteps[timestep])}")
- print("electrons_across_timesteps:", electrons_across_timesteps)
-
- _results.append(_result)
- print(f"Completed timestep {timestep}")
-
- except Exception as e:
- print(f"Error processing timestep {timestep}: {str(e)}")
- raise
-
- # Final debug output
- print("\nFinal electrons_across_timesteps:")
- print(f"Shape: {electrons_across_timesteps.shape}")
- print(f"Min value: {np.min(electrons_across_timesteps)}")
- print(f"Max value: {np.max(electrons_across_timesteps)}")
- print(f"Mean value: {np.mean(electrons_across_timesteps)}")
-
- # Plot results
- plot_detector_results1(electrons_across_timesteps)
-
- return _results
-
-def plot_detector_results(electrons_per_detector_time):
- """Plot electron counts for each detector across timesteps"""
- plt.figure(figsize=(12, 8))
- timesteps = range(len(electrons_per_detector_time))
-
- # Plot for each detector
- for detector in range(4):
- plt.plot(timesteps, electrons_per_detector_time[:, detector],
- label=f'Detector {detector+1}',
- marker='o')
-
- plt.title("Electrons Captured by Detectors Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Number of Electrons")
- plt.legend()
- plt.grid(True)
-
- # Add value annotations
- for detector in range(4):
- for t, value in zip(timesteps, electrons_per_detector_time[:, detector]):
- plt.annotate(f'{value:.2e}',
- (t, value),
- textcoords="offset points",
- xytext=(0,10),
- ha='center',
- rotation=45)
-
- plt.tight_layout()
- plt.savefig("detector_electrons_over_time(VDB-multi).png", dpi=300, bbox_inches='tight')
- plt.close()
-
- # Create additional visualization - normalized values
- plt.figure(figsize=(12, 8))
- normalized_electrons = electrons_per_detector_time / np.max(electrons_per_detector_time)
-
- for detector in range(4):
- plt.plot(timesteps, normalized_electrons[:, detector],
- label=f'Detector {detector+1}',
- marker='o')
-
- plt.title("Normalized Electrons Captured by Detectors Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Normalized Electron Count")
- plt.legend()
- plt.grid(True)
- plt.tight_layout()
- plt.savefig("detector_electrons_over_time_normalized(VDB-multi-000003).png", dpi=300, bbox_inches='tight')
- plt.close()
-
-def plot_detector_results1(electrons_per_detector_time):
- """
- Plot electron counts for detector pairs across timesteps
-
- Args:
- electrons_per_detector_time: Array of shape (num_timesteps, 4) containing
- electron counts for each detector at each timestep
- :param electrons_per_detector_time:
- """
- # Create figure for detectors 1 and 2
- plt.figure(figsize=(12, 8))
- timesteps = range(len(electrons_per_detector_time))
-
- # Plot detectors 1 and 2
- plt.plot(timesteps, electrons_per_detector_time[:, 0],
- label='Detector 1', marker='o', linewidth=1, markevery=20)
- plt.plot(timesteps, electrons_per_detector_time[:, 1],
- label='Detector 2', marker='s', linewidth=1, markevery=20)
-
- plt.title("Electrons Captured by Detectors 1 and 2 Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Number of Electrons")
- plt.legend()
- plt.grid(True)
- plt.tight_layout()
- plt.savefig("detectors_1_2_over_time(VDB-multi-000003).png", dpi=300, bbox_inches='tight')
- plt.close()
-
- # Create figure for detectors 3 and 4
- plt.figure(figsize=(12, 8))
-
- # Plot detectors 3 and 4
- plt.plot(timesteps, electrons_per_detector_time[:, 2],
- label='Detector 3', marker='o', linewidth=1, markevery=20)
- plt.plot(timesteps, electrons_per_detector_time[:, 3],
- label='Detector 4', marker='s', linewidth=1, markevery=20)
-
- plt.title("Electrons Captured by Detectors 3 and 4 Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Number of Electrons")
- plt.legend()
- plt.grid(True)
- plt.tight_layout()
- plt.savefig("detectors_3_4_over_time(VDB-multi-000003).png", dpi=300, bbox_inches='tight')
- plt.close()
-
-if __name__ == "__main__":
- results = main()
\ No newline at end of file
diff --git a/apps/showcases/FreeSurface/raytracerVDB-multi.py b/apps/showcases/FreeSurface/raytracerVDB-multi.py
deleted file mode 100644
index 6c807e4f8011310775c9d5baa94e3fd3cbb4eb83..0000000000000000000000000000000000000000
--- a/apps/showcases/FreeSurface/raytracerVDB-multi.py
+++ /dev/null
@@ -1,1798 +0,0 @@
-import pyopenvdb as vdb
-import os
-from math import floor
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib.colors as colors
-from matplotlib.animation import FuncAnimation
-from typing import Tuple, List
-import scipy.ndimage as ndimage
-from scipy.ndimage import gaussian_filter
-from scipy.spatial.transform import Rotation as Ro
-from scipy.spatial.transform import Slerp
-from scipy.interpolate import PchipInterpolator
-from sympy.printing.tree import print_node
-
-
-def runRayTracer(x_num, y_num, z_num, fill_levels, current_timestep):
-
- print("current_timestep", current_timestep)
- # fill_levels = np.clip(fill_levels, 0, 1)
- fill_levels = fill_levels.astype(np.float64)
-
- # Define the tolerance
- tolerance = 1e-8
-
- # Check if any value in fill_levels is outside the [0, 1] range
- out_of_range = np.logical_or(fill_levels < -tolerance, fill_levels > 1 + tolerance)
- if np.any(out_of_range):
- out_of_range_values = fill_levels[out_of_range]
- out_of_range_indices = np.where(out_of_range)
-
- for idx, value in zip(zip(*out_of_range_indices), out_of_range_values):
- print(f"Index: {idx}, Value: {value}")
-
- # raise ValueError(f"Fill level values must be between 0 and 1 inclusive (with tolerance of {tolerance}).")
-
- # Check for fill levels within the range [0, 1]
- in_range = np.logical_and(fill_levels > 0, fill_levels < 1)
- in_range_values = fill_levels[in_range]
- in_range_indices = np.where(in_range)
-
- # print("Fill level values within range (0 to 1):")
- # for idx, value in zip(zip(*in_range_indices), in_range_values):
- # if idx == (np.int64(309),):
- # print(f"Index: {idx}, Value: {value}")
-
- # Constants
- pi: float = np.pi # Mathematical constant π
- e: float = 1.60217662e-19 # Elementary charge in Coulombs
- beam_current: float = 3e-3 # Beam current in Amperes
- dwell_time: float = 0.32e-6 # Dwell time in seconds
- A: float = 21.5 # Approximation for Ti6Al4V
- # Ti6Al4V = Ti6Al4V.create_Ti6Al4V(T)
- # A: float = Ti6Al4V.atomic_number
- # print('Ti6Al4V.elements: ', Ti6Al4V.elements)
- # print('Ti6Al4V.composition: ', Ti6Al4V.composition)
- # print('A: ', A)
- dx: float = 3e-6 # Cell size in meters
- dy: float = dx # Cell size in meters
- dz: float = dx # Cell size in meters
- threshold: float = 0 # Beam profile threshold
- # The voltage at detector d can be computed using the transimpedance T (in units of V /A ) of a virtual amplifier in the context of the model
- transimpedance: float = 1e5 # Transimpedance in Volts per Ampere
- eV_threshold = 50 # Threshold energy in eV
- threshold_energy = eV_threshold * e # Threshold energy in joules
- # Define the work distance (height of the electron beam gun)
- # w_d = 12e-3 # Work distance in meters (272 mm)
-
- # Domain setup (Section 2.1.1)
- # Define the 3D domain for the simulation
- x_min, x_max = -0.5 * x_num * dx, 0.5 * x_num * dx # -10 mm to 10 mm
- # print('x_min, x_max:', x_min, x_max)
- y_min, y_max = -0.5 * y_num * dy, 0.5 * y_num * dy # -10 mm to 10 mm
- # print('y_min, y_max:', y_min, y_max)
- z_min, z_max = 0, z_num * dz
- # print('z_min, z_max:', z_min, z_max)
- # x_num = int((x_max - x_min) / dx) + 1
- # y_num = int((y_max - y_min) / dy) + 1
- # z_num = int((z_max - z_min) / dz) + 1
- # print(f"Domain dimensions: x={x_num}, y={y_num}, z={z_num}")
- # domain = np.zeros((z_num, y_num, x_num)) # (z, y , x)
- # print('domain shape: ', np.shape(domain))
- # domain[0:2, :, :] = 1 # Solid cells
- # domain[2, :, :] = 0.5 # Interface layer
- fill_levels_reshaped = np.reshape(fill_levels, (z_num, y_num, x_num))
- # fill_levels_reshaped = fill_levels
- # print('fill_levels_reshaped:', fill_levels_reshaped)
- out_of_range_fill_levels = np.logical_or(fill_levels_reshaped < -tolerance, fill_levels_reshaped > 1 + tolerance)
- if np.any(out_of_range_fill_levels):
- out_of_range_fill_level_values = fill_levels_reshaped[out_of_range_fill_levels]
- out_of_range_fill_level_indices = np.where(out_of_range_fill_levels)
-
- for idx, value in zip(zip(*out_of_range_fill_level_indices), out_of_range_fill_level_values):
- print(f"Index: {idx}, Value: {value}")
- raise ValueError(f"Fill level values must be between 0 and 1 inclusive (with tolerance of {tolerance}).")
-
- sigma = 0.0
- # Different sigma for each dimension
- # sigma = [0.0, 1.0, 0.5] # [z, y, x]
- fill_levels_smoothed = ndimage.gaussian_filter(fill_levels_reshaped, sigma=sigma, mode='nearest', radius=1)
-
- fill_levels_smoothed = np.clip(
- fill_levels_smoothed,
- a_min=np.float64(0.0),
- a_max=np.float64(1.0)
- )
-
- # Check if values are still in bounds
- if not (np.all(fill_levels_smoothed >= 0) and np.all(fill_levels_smoothed <= 1)):
- raise ValueError("Smoothed fill levels out of bounds!")
-
-
- def interpolate_normals(_normal_below: np.ndarray, _normal_above: np.ndarray, t: float) -> np.ndarray:
- """
- Interpolate between two normals using spherical linear interpolation (SLERP)
- """
- # Input validation with better error messages
- if not (0 <= t <= 1):
- raise ValueError(f"Interpolation parameter t must be in [0,1], got {t}")
- if not (isinstance(_normal_below, np.ndarray) and isinstance(_normal_above, np.ndarray)):
- raise ValueError("Normals must be numpy arrays")
-
- # Ensure normals are normalized and correct type
- _normal_below = normalize_vector(_normal_below.astype(np.float64))
- _normal_above = normalize_vector(_normal_above.astype(np.float64))
-
- # Compute dot product with numerical stability
- dot_product = np.clip(np.dot(_normal_below, _normal_above), -1.0, 1.0)
-
- # Better thresholds
- PARALLEL_THRESHOLD = 0.9999 # cos(1°) ≈ 0.9998
- ANTIPARALLEL_THRESHOLD = -0.9999
-
- # Case 1: Normals nearly parallel
- if dot_product > PARALLEL_THRESHOLD:
- _normal_interp = (1.0 - t) * _normal_below + t * _normal_above
- return normalize_vector(_normal_interp)
-
- # Case 2: Normals nearly antiparallel
- if dot_product < ANTIPARALLEL_THRESHOLD:
- # Find perpendicular vector more robustly
- perp = np.array([1.0, 0.0, 0.0])
- if abs(np.dot(_normal_below, perp)) > 0.9:
- perp = np.array([0.0, 1.0, 0.0])
- if abs(np.dot(_normal_below, perp)) > 0.9:
- perp = np.array([0.0, 0.0, 1.0])
- perp = normalize_vector(np.cross(_normal_below, perp))
-
- # Rotate 180 degrees around perp
- angle = np.pi * t
- _normal_interp = _normal_below * np.cos(angle) + perp * np.sin(angle)
- return normalize_vector(_normal_interp)
-
- # Case 3: Standard SLERP with better numerical stability
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
-
- if abs(sin_theta) < 1e-10:
- return normalize_vector(_normal_below)
-
- scale0 = np.sin((1.0 - t) * theta) / sin_theta
- scale1 = np.sin(t * theta) / sin_theta
- _normal_interp = scale0 * _normal_below + scale1 * _normal_above
-
- return normalize_vector(_normal_interp)
-
-
- def interpolate_fill_level_and_normal(_fill_level, surface_normals, _y, _x):
- # Validate input
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
- # Find the indices where the fill level crosses 0.5 in either direction
- indices_above = np.where((_fill_level[:, _y, _x][:-1] > 0.5) & (_fill_level[:, _y, _x][1:] <= 0.5))[0]
- print('indices_above:', indices_above)
- indices_below = np.where((_fill_level[:, _y, _x][:-1] <= 0.5) & (_fill_level[:, _y, _x][1:] > 0.5))[0]
- print('indices_below:', indices_below)
-
- # if len(indices_above) == 0 and len(indices_below) == 0:
- # Fill level does not cross 0.5 at the given (_x, _y) coordinate
- # return None, None, None, None, None, None, None
-
- if len(indices_above) > 0:
- print("len(indices_above) > 0")
- # Fill level crosses 0.5 from above to below
- z_index = indices_above[0]
- print('z_index:', z_index)
- _fill_level_above = _fill_level[z_index, _y, _x]
- print('_fill_level_above:', _fill_level_above)
- _fill_level_below = _fill_level[z_index + 1, _y, _x]
- print('_fill_level_below:', _fill_level_below)
- elif len(indices_below) > 0:
- print("len(indices_below) > 0")
- # Fill level crosses 0.5 from below to above
- z_index = indices_below[0]
- print('z_index:', z_index)
- _fill_level_below = _fill_level[z_index, _y, _x]
- print('_fill_level_below:', _fill_level_below)
- _fill_level_above = _fill_level[z_index + 1, _y, _x]
- print('_fill_level_above:', _fill_level_above)
- else:
- # Fill level does not cross 0.5 at the given (_x, _y) coordinate
- print("No interface cells found")
- return None, None, None, None, None, None, None
-
- # Interpolate the z-coordinate where the fill level is 0.5
- denominator = _fill_level_below - _fill_level_above
- if abs(denominator) < 1e-6:
- # Handle the case where the denominator is very close to zero
- _z_interp = z_index + 0.5
- else:
- _z_interp = z_index + ((0.5 - _fill_level_above) / denominator)
- # _z_interp = z_index + (0.5 - bilinear_interpolate(_x, _y, z_index, _fill_level)) / (bilinear_interpolate(_x, _y, z_index+1, _fill_level) - bilinear_interpolate(_x, _y, z_index, _fill_level))
- print('_z_interp:', _z_interp)
-
- _z_height = (_z_interp + 0.5) * dz
- print('_z_height:', _z_height)
-
- # Interpolate the surface normal at the interpolated z-coordinate
- _normal_below = surface_normals[z_index, _y, _x].astype(np.float64)
- print('_normal_below:', _normal_below)
- _normal_above = surface_normals[z_index + 1, _y, _x].astype(np.float64)
- print('_normal_above:', _normal_above)
-
- t = _z_interp - z_index
-
- # Use spherical interpolation (slerp) for more accurate normal interpolation
- dot_product = np.dot(_normal_below, _normal_above).astype(np.float64)
- if dot_product > np.cos(np.deg2rad(1)): # Threshold for using linear interpolation
- # If normals are very close, use linear interpolation
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- # Use spherical interpolation
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + (np.sin(t*theta) / sin_theta) * _normal_above
- # _normal_interp = quaternion_interpolate(_normal_below, _normal_above, t)
-
- _normal_interp = normalize_vector(_normal_interp)
- print('_normal_interp:', _normal_interp)
-
- return _z_interp, _z_height, _fill_level_below, _fill_level_above, _normal_interp, _normal_below, _normal_above
-
-
- def interpolate_fill_level_and_normal1(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- # Validate input
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at the given (x,y) position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom (reverse order)
- for z in range(len(column)-2, -1, -1): # Start from second-to-last index
- fill_current = column[z]
- fill_next = column[z + 1]
-
- # Check if we cross 0.5 threshold
- if (fill_current - 0.5) * (fill_next - 0.5) <= 0:
- # Found interface
- z_index = z
- _fill_level_above = fill_current
- _fill_level_below = fill_next
-
- print(f"Found interface at z={z_index}")
- print(f"fill_level_above: {_fill_level_above}")
- print(f"fill_level_below: {_fill_level_below}")
-
- # Interpolate z position
- denominator = _fill_level_below - _fill_level_above
- if abs(denominator) < 1e-6:
- _z_interp = z_index + 0.5
- else:
- _z_interp = z_index + ((0.5 - _fill_level_above) / denominator)
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and interpolate normals
- _normal_below = surface_normals[z_index, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z_index + 1, _y, _x].astype(np.float64)
-
- # Interpolate normals
- t = _z_interp - z_index
- dot_product = np.dot(_normal_below, _normal_above)
-
- if dot_product > np.cos(np.deg2rad(1)):
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + \
- (np.sin(t*theta) / sin_theta) * _normal_above
-
- _normal_interp = normalize_vector(_normal_interp)
-
- return _z_interp, _z_height, _fill_level_below, _fill_level_above, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
-
- def interpolate_fill_level_and_normal2(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at the given (x,y) position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom
- for z in range(len(column)-2, -1, -1):
- fill_current = column[z]
- fill_next = column[z + 1]
-
- # Look for significant transition from liquid (>0.5) to gas (<0.2)
- if fill_current > 0.5 and fill_next < 0.2:
- print(f"Found interface at z={z}")
- print(f"fill_level_above: {fill_next}") # gas phase
- print(f"fill_level_below: {fill_current}") # liquid phase
-
- # Interpolate z position
- denominator = fill_next - fill_current
- if abs(denominator) < 1e-6:
- _z_interp = z + 0.5
- else:
- _z_interp = z + ((0.5 - fill_current) / denominator)
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and interpolate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- # Interpolate normal
- t = _z_interp - z
- dot_product = np.dot(_normal_below, _normal_above)
-
- if dot_product > np.cos(np.deg2rad(1)):
- _normal_interp = _normal_below + t * (_normal_above - _normal_below)
- else:
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
- _normal_interp = (np.sin((1-t)*theta) / sin_theta) * _normal_below + \
- (np.sin(t*theta) / sin_theta) * _normal_above
-
- _normal_interp = normalize_vector(_normal_interp)
-
- return _z_interp, _z_height, fill_next, fill_current, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
-
- def interpolate_fill_level_and_normal3(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals, starting from top"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals, {surface_normals.ndim}")
-
- # Get fill levels at position
- column = _fill_level[:, _y, _x]
-
- # Scan from top to bottom
- for z in range(len(column)-2, -1, -1):
- current_fill = column[z]
- next_fill = column[z + 1]
-
- # Look for liquid-gas interface (current > 0, next = 0)
- if current_fill > 0 and next_fill == 0:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}") # liquid phase
- print(f"fill_level_above: {next_fill}") # gas phase
-
- # Interface is at z since this is the last liquid cell
- denominator = next_fill - current_fill
- if abs(denominator) < 1e-6:
- _z_interp = z + 0.5
- else:
- _z_interp = z + ((0.5 - current_fill) / denominator)
-
- print('_z_interp', _z_interp)
- _z_height = (_z_interp + 0.5) * dz
- print("_z_height:", _z_height)
-
- # Get and interpolate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- # Interpolate normal
- t = _z_interp - z
-
- # Use default normal if either normal is invalid
- if np.any(np.isnan(_normal_below)) or np.any(np.isnan(_normal_above)):
- _normal_interp = np.array([0.0, 0.0, 1.0])
- else:
- # Interpolate normal at interface
- _normal_interp = interpolate_normals(_normal_below, _normal_above, t)
-
- # _normal_interp = interpolate_normals(_normal_below, _normal_above, t)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
-
- def interpolate_fill_level_and_normal4(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals: {surface_normals.ndim}")
-
- # Get fill levels at position
- column = _fill_level[:, _y, _x]
-
- # Constants
- INTERFACE_THRESHOLD = 0.5 # Standard interface threshold
- GRADIENT_THRESHOLD = 0.1 # Minimum gradient for interface detection
- TOLERANCE = 1e-6 # Numerical tolerance
-
- # Scan from top to bottom
- for z in range(len(column)-2, -1, -1):
- current_fill = column[z]
- next_fill = column[z + 1]
-
- # Check for significant gradient and crossing of 0.5
- gradient = abs(next_fill - current_fill)
- if gradient > GRADIENT_THRESHOLD:
- # Check if we cross the interface threshold
- if ((current_fill > INTERFACE_THRESHOLD > next_fill) or
- (current_fill < INTERFACE_THRESHOLD < next_fill)):
-
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}")
- print(f"fill_level_above: {next_fill}")
-
- # Calculate interpolation parameter
- denominator = next_fill - current_fill
- if abs(denominator) < TOLERANCE:
- t_interface = 0.5
- _z_interp = z + 0.5
- else:
- t_interface = (INTERFACE_THRESHOLD - current_fill) / denominator
- t_interface = np.clip(t_interface, 0.0, 1.0) # Ensure t is in [0,1]
- _z_interp = z + t_interface
-
- print('_z_interp', _z_interp)
- _z_height = (_z_interp + 0.5) * dz
-
- # Get normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- # Handle invalid normals
- if np.any(np.isnan(_normal_below)) or np.any(np.isnan(_normal_above)):
- _normal_interp = np.array([0.0, 0.0, 1.0])
- else:
- # Use clipped t_interface for normal interpolation
- _normal_interp = interpolate_normals(_normal_below, _normal_above, t_interface)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- print("No interface found")
- return None, None, None, None, None, None, None
-
- def interpolate_fill_level_and_normal5(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals: {surface_normals.ndim}")
-
- # Get fill levels at position
- _column = _fill_level[:, _y, _x]
-
- # Constants with physical meaning
- INTERFACE_THRESHOLD = 0.5 # Standard VOF interface threshold
- GRADIENT_THRESHOLD = 0.1 # Minimum gradient for interface detection
- TOLERANCE = 1e-6 # Numerical tolerance
-
- # Scan from top to bottom
- for z in range(len(_column)-2, -1, -1):
- current_fill = _column[z]
- next_fill = _column[z + 1]
-
- # Check for significant gradient
- gradient = abs(next_fill - current_fill)
- if gradient > GRADIENT_THRESHOLD:
- # Check for interface crossing (VOF method)
- crosses_interface = (
- (current_fill >= INTERFACE_THRESHOLD >= next_fill) or
- (current_fill <= INTERFACE_THRESHOLD <= next_fill)
- )
-
- if crosses_interface:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}")
- print(f"fill_level_above: {next_fill}")
- # Linear interpolation with better numerical stability
- denominator = next_fill - current_fill
- if abs(denominator) < TOLERANCE:
- _z_interp = z + 0.5
- t_interface = 0.5
- else:
- # Calculate interpolation parameter
- t_interface = (INTERFACE_THRESHOLD - current_fill) / denominator
- if not (0 <= t_interface <= 1):
- error_msg = (
- f"\nInterpolation parameter t_interface={t_interface:.6f} is out of bounds [0,1]\n"
- f"Debug information:\n"
- f" - z-index: {z}\n"
- f" - current_fill: {current_fill:.6f}\n"
- f" - next_fill: {next_fill:.6f}\n"
- f" - denominator: {denominator:.6f}\n"
- f" - INTERFACE_THRESHOLD: {INTERFACE_THRESHOLD}\n"
- f" - Position (y,x): ({_y}, {_x})\n"
- f" - Fill level values around interface:\n"
- )
- # Nearby fill levels for context
- for i in range(max(0, z-2), min(len(_column), z+3)):
- error_msg += f" {i:3d} | {_column[i]:.6f}\n"
-
- raise ValueError(error_msg)
- # t_interface = np.clip(t_interface, 0.0, 1.0)
- _z_interp = z + t_interface
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and validate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- _normal_interp = interpolate_normals(_normal_below, _normal_above, t_interface)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- return None, None, None, None, None, None, None
-
-
- def compute_surface_normals_sobel(_fill_levels: np.ndarray) -> np.ndarray:
- """
- Compute surface normals using Sobel filter with enhanced numerical precision.
- """
- # Convert input to high precision at the start
- _fill_levels = _fill_levels.astype(np.float64)
-
- # Input validation
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- # Define tolerances for different checks
- TOLERANCE = {
- 'fill_level_bounds': 1e-10,
- 'zero': 1e-10,
- 'norm': 1e-10,
- 'normalization': 1e-8
- }
-
- # Validate domain bounds with high precision
- '''invalid_elements = np.logical_or(_fill_levels < -TOLERANCE['fill_level_bounds'],
- _fill_levels > 1 + TOLERANCE['fill_level_bounds'])
- if np.any(invalid_elements):
- invalid_indices = np.where(invalid_elements)
- for _z, _y, _x in zip(*invalid_indices):
- _value = _fill_levels[_z, _y, _x]
- print(f"Warning: Invalid value {_value} at [z={_z}, y={_y}, x={_x}]")'''
-
- # Compute gradients using Sobel filter with explicit high precision
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64)
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64)
- _grad_z = -1 * ndimage.sobel(_fill_levels, axis=0, output=np.float64) # multiply by -1 to invert the surface normals to point in +z direction
-
- # Stack gradients into normal vector array with maintained precision
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Compute norm with high precision
- norm = np.linalg.norm(_n_S, axis=-1)
-
- # Create masks for valid regions
- interface_mask = (_fill_levels > TOLERANCE['zero']) & (_fill_levels < 1 - TOLERANCE['zero'])
- nonzero_mask = norm > TOLERANCE['norm']
- valid_mask = interface_mask & nonzero_mask
-
- # Normalize vectors only where valid
- # _n_S[valid_mask] /= norm[valid_mask][..., np.newaxis]
- # _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized = np.full_like(_n_S, np.nan)
- _n_S_normalized[valid_mask] = _n_S[valid_mask] / norm[valid_mask, np.newaxis]
-
- # Verify normalization
- if not np.allclose(np.linalg.norm(_n_S_normalized[valid_mask], axis=-1),1.0, atol=TOLERANCE['normalization']):
- raise ValueError("Normalization failed for some vectors")
-
- # return _n_S
- return _n_S_normalized
-
-
- def compute_surface_normals_sobel1(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'fill_level_bounds': 1e-10,
- 'zero': 1e-10,
- 'norm': 1e-10,
- 'normalization': 1e-8
- }
-
- # Compute gradients
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64)
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64)
- _grad_z = -1 * ndimage.sobel(_fill_levels, axis=0, output=np.float64)
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Compute norm
- norm = np.linalg.norm(_n_S, axis=-1)
-
- # Create masks for valid regions:
- # 1. Interface cells (fill level between 0 and 1)
- interface_mask = (_fill_levels > TOLERANCE['zero']) & (_fill_levels < 1 - TOLERANCE['zero'])
-
- # 2. Cells adjacent to fill level transitions
- # filled_mask = _fill_levels >= 1 - TOLERANCE['zero']
- filled_mask = _fill_levels >= 1
- # empty_mask = _fill_levels <= TOLERANCE['zero']
- empty_mask = _fill_levels <= 0
-
- # Check for transitions in z direction
- z_transitions = np.zeros_like(_fill_levels, dtype=bool)
- z_transitions[:-1, :, :] |= (filled_mask[:-1, :, :] & empty_mask[1:, :, :])
- z_transitions[1:, :, :] |= (empty_mask[:-1, :, :] & filled_mask[1:, :, :])
-
- # Combine masks
- valid_mask = interface_mask | z_transitions
-
- # Default normal for sharp transitions (pointing up)
- _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized[..., 2] = 1.0 # Default z-direction normal
-
- # Expand valid_mask to match _n_S shape
- valid_mask_expanded = valid_mask[..., np.newaxis]
-
- # Where we have valid gradients, use them
- gradient_mask = (norm > TOLERANCE['norm'])[..., np.newaxis]
- mask = valid_mask_expanded & gradient_mask
-
- # Normalize vectors where mask is True
- norm_expanded = norm[..., np.newaxis]
- _n_S_normalized = np.where(mask, _n_S / np.where(norm_expanded > 0, norm_expanded, 1.0), _n_S_normalized)
-
- return _n_S_normalized
-
-
- def compute_surface_normals_sobel2(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'fill_level_bounds': 1e-10,
- 'zero': 1e-10,
- 'norm': 1e-10,
- 'normalization': 1e-8
- }
-
- # Compute gradients
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64)
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64)
- _grad_z = -1 * ndimage.sobel(_fill_levels, axis=0, output=np.float64)
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Compute norm
- norm = np.linalg.norm(_n_S, axis=-1)
-
- # Create masks for valid regions:
- # 1. Interface cells (fill level between 0 and 1)
- interface_mask = (_fill_levels > TOLERANCE['zero']) & (_fill_levels < 1 - TOLERANCE['zero'])
-
- # 2. Cells adjacent to fill level transitions
- filled_mask = _fill_levels >= 1 - TOLERANCE['zero']
- empty_mask = _fill_levels <= TOLERANCE['zero']
-
- # Check for transitions in z direction
- z_transitions = np.zeros_like(_fill_levels, dtype=bool)
- z_transitions[:-1, :, :] |= (filled_mask[:-1, :, :] & empty_mask[1:, :, :])
- z_transitions[1:, :, :] |= (empty_mask[:-1, :, :] & filled_mask[1:, :, :])
-
- # Combine masks
- valid_mask = interface_mask | z_transitions
-
- # Initialize output array with NaN
- _n_S_normalized = np.full_like(_n_S, np.nan)
-
- # Create mask for valid gradients
- gradient_mask = norm > TOLERANCE['norm']
-
- # Combine masks and reshape for broadcasting
- mask = valid_mask & gradient_mask
-
- # Normalize only where mask is True
- for i in range(3): # For each component (x, y, z)
- _n_S_normalized[..., i] = np.where(mask,
- _n_S[..., i] / np.where(norm > 0, norm, 1.0),
- np.nan)
-
- return _n_S_normalized
-
- def compute_surface_normals_sobel3(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with default upward normal
- _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized[..., 2] = 1.0
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
-
- def interpolate_normals_new(_normal_below: np.ndarray, _normal_above: np.ndarray, t: float) -> np.ndarray:
- """Interpolate between two normals using SLERP"""
- if not (0 <= t <= 1):
- raise ValueError(f"Interpolation parameter t must be in [0,1], got {t}")
-
- # Ensure normalized vectors
- _normal_below = normalize_vector(_normal_below.astype(np.float64))
- _normal_above = normalize_vector(_normal_above.astype(np.float64))
-
- # Compute dot product
- dot_product = np.clip(np.dot(_normal_below, _normal_above), -1.0, 1.0)
-
- # Thresholds for special cases
- PARALLEL_THRESHOLD = 0.9999
- ANTIPARALLEL_THRESHOLD = -0.9999
-
- # Case 1: Nearly parallel
- if dot_product > PARALLEL_THRESHOLD:
- return normalize_vector((1.0 - t) * _normal_below + t * _normal_above)
-
- # Case 2: Nearly antiparallel
- if dot_product < ANTIPARALLEL_THRESHOLD:
- perp = np.array([1.0, 0.0, 0.0])
- if abs(np.dot(_normal_below, perp)) > 0.9:
- perp = np.array([0.0, 1.0, 0.0])
- perp = normalize_vector(np.cross(_normal_below, perp))
- angle = np.pi * t
- return normalize_vector(_normal_below * np.cos(angle) + perp * np.sin(angle))
-
- # Case 3: Standard SLERP
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
-
- if abs(sin_theta) < 1e-10:
- return normalize_vector(_normal_below)
-
- scale0 = np.sin((1.0 - t) * theta) / sin_theta
- scale1 = np.sin(t * theta) / sin_theta
- return normalize_vector(scale0 * _normal_below + scale1 * _normal_above)
-
-
- def interpolate_fill_level_and_normal_new(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals: {surface_normals.ndim}")
-
- # Get fill levels at position
- _column = _fill_level[:, _y, _x]
-
- # Constants with physical meaning
- INTERFACE_THRESHOLD = 0.5 # Standard VOF interface threshold
- GRADIENT_THRESHOLD = 0.1 # Minimum gradient for interface detection
- TOLERANCE = 1e-6 # Numerical tolerance
-
- # Scan from top to bottom
- for z in range(len(_column)-2, -1, -1):
- current_fill = _column[z]
- next_fill = _column[z + 1]
-
- # Check for significant gradient
- gradient = abs(next_fill - current_fill)
- if gradient > GRADIENT_THRESHOLD:
- # Check for interface crossing (VOF method)
- '''crosses_interface = (
- (current_fill >= INTERFACE_THRESHOLD - TOLERANCE >= next_fill) or
- (current_fill <= INTERFACE_THRESHOLD + TOLERANCE <= next_fill)
- )'''
- crosses_interface = (
- (current_fill >= (INTERFACE_THRESHOLD - TOLERANCE) and
- next_fill <= (INTERFACE_THRESHOLD + TOLERANCE))
- )
-
- if crosses_interface:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}")
- print(f"fill_level_above: {next_fill}")
- # Linear interpolation with better numerical stability
- denominator = next_fill - current_fill
- if abs(denominator) < TOLERANCE:
- _z_interp = z + 0.5
- t_interface = 0.5
- else:
- # Calculate interpolation parameter
- t_interface = (INTERFACE_THRESHOLD - current_fill) / denominator
- if not (0 <= t_interface <= 1):
- error_msg = (
- f"\nInterpolation parameter t_interface={t_interface:.6f} is out of bounds [0,1]\n"
- f"Debug information:\n"
- f" - z-index: {z}\n"
- f" - current_fill: {current_fill:.6f}\n"
- f" - next_fill: {next_fill:.6f}\n"
- f" - denominator: {denominator:.6f}\n"
- f" - INTERFACE_THRESHOLD: {INTERFACE_THRESHOLD}\n"
- f" - Position (y,x): ({_y}, {_x})\n"
- f" - Fill level values around interface:\n"
- )
- # Nearby fill levels for context
- for i in range(max(0, z-2), min(len(_column), z+3)):
- error_msg += f" {i:3d} | {_column[i]:.6f}\n"
-
- raise ValueError(error_msg)
- # t_interface = np.clip(t_interface, 0.0, 1.0)
- _z_interp = z + t_interface
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and validate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- _normal_interp = interpolate_normals_new(_normal_below, _normal_above, t_interface)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- return None, None, None, None, None, None, None
-
-
- def compute_surface_normals_sobel_zeros(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with default upward normal
- _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized[..., 2] = 1.0
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
-
- def compute_surface_normals_sobel_nans(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter with NaN initialization"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with NaN
- _n_S_normalized = np.full_like(_n_S, np.nan)
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
-
- # Compute surface normals (Section 2.1.1)
- # n_S, (grad_x, grad_y, grad_z) = compute_surface_normals(domain)
- n_S = compute_surface_normals_sobel_nans(fill_levels_smoothed)
- # print('n_S: ', np.shape(n_S))
-
- # Beam profile (Section 2.1)
- # Define the Gaussian beam profile
- sigma: float = 500e-6 / (4 * dx) # Beam profile sigma
- x = np.arange(-np.floor(min(x_num//4, 2*sigma)), 1+np.floor(min((x_num//4), 2*sigma)))
- # print('x:', x)
- y = np.arange(-floor(min(y_num//4, 2*sigma)), 1+floor(min((y_num//4), 2*sigma)))
- # print('y:', y)
- X, Y = np.meshgrid(x, y)
- beam_profile = np.exp(-(X ** 2 + Y ** 2) / (2 * sigma ** 2))
- beam_profile /= np.sum(beam_profile) # Normalize beam profile
- print('beam_profile shape: ', np.shape(beam_profile))
-
- beam_size = np.shape(beam_profile)[0] * np.shape(beam_profile)[1]
- # print('Beam size:', beam_size)
-
- # Detector setup (Section 2.1.1)
- # Define positions of the detectors
- detector_height = 4 * z_max # Detector height in meters (272 mm)
- detector_spacing_x = 1.5 * dx * x_num # Detector spacing in meters (127 mm)
- detector_spacing_y = detector_spacing_x # 1.5 * dy * y_num # Detector spacing in meters (127 mm)
- detector_diameter = 4 * dx # Detector diameter in meters (50 mm)
-
- # Define detector positions relative to the origin
- detector_positions = np.array([
- [0.5 * detector_spacing_x, 0, detector_height],
- [-0.5 * detector_spacing_x, 0, detector_height],
- [0, 0.5 * detector_spacing_y, detector_height],
- [0, -0.5 * detector_spacing_y, detector_height]
- ], dtype=np.float64)
- # print('detector_positions:', detector_positions)
-
-
- # Define detector vertices for solid angle calculation
- def get_detector_vertices(detector_position: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
- # Calculate the side length of the detector
- side_length : float = np.sqrt(3) * (detector_diameter / 2)
- height : float = np.sqrt(3) * (side_length / 2)
-
- # Determine the orientation based on the detector position
- if np.allclose(detector_position[:2], [0, 0]):
- raise ValueError("Detector position cannot be at the origin (0, 0, 0)")
- elif detector_position[0] > 0: # Right detector
- _v1 = detector_position + np.array([2 / 3 * height, 0, 0], dtype=np.float64)
- _v2 = detector_position + np.array([-height / 3, side_length / 2, 0], dtype=np.float64)
- _v3 = detector_position + np.array([-height / 3, -side_length / 2, 0], dtype=np.float64)
- elif detector_position[0] < 0: # Left detector
- _v1 = detector_position + np.array([-2 / 3 * height, 0, 0], dtype=np.float64)
- _v2 = detector_position + np.array([height / 3, -side_length / 2, 0], dtype=np.float64)
- _v3 = detector_position + np.array([height / 3, side_length / 2, 0], dtype=np.float64)
- elif detector_position[1] > 0: # Back detector
- _v1 = detector_position + np.array([0, 2 / 3 * height, 0], dtype=np.float64)
- _v2 = detector_position + np.array([-side_length / 2, -height / 3, 0], dtype=np.float64)
- _v3 = detector_position + np.array([side_length / 2, -height / 3, 0], dtype=np.float64)
- elif detector_position[1] < 0: # Front detector
- _v1 = detector_position + np.array([0, -2 / 3 * height, 0], dtype=np.float64)
- _v2 = detector_position + np.array([side_length / 2, height / 3, 0], dtype=np.float64)
- _v3 = detector_position + np.array([-side_length / 2, height / 3, 0], dtype=np.float64)
- else:
- raise ValueError("Invalid detector position. Must be in one of the four expected positions.")
-
- return _v1, _v2, _v3
-
-
- # Compute solid angle using Equation A1 (Appendix A)
- def compute_solid_angle(vector1: np.ndarray, vector2: np.ndarray, vector3: np.ndarray) -> float:
- # Input validation
- if not isinstance(vector1, np.ndarray) or not isinstance(vector2, np.ndarray) or not isinstance(vector3, np.ndarray):
- raise ValueError("Inputs must be numpy arrays")
- if vector1.shape != (3,) or vector2.shape != (3,) or vector3.shape != (3,):
- raise ValueError("Inputs must be 3D vectors")
- if not np.isfinite(vector1).all() or not np.isfinite(vector2).all() or not np.isfinite(vector3).all():
- raise ValueError("Inputs must be finite")
-
- # Convert to high precision once
- vectors = [vec.astype(np.float64) for vec in (vector1, vector2, vector3)]
- cross_product = np.cross(vectors[1], vectors[2])
- # Check if the cross product is nearly zero
- if np.linalg.norm(cross_product) < 1e-16:
- raise ValueError("Vectors are nearly coplanar")
-
- # Compute solid angle
- try:
- # Compute the numerator of the solid angle formula
- numerator: float = np.dot(vectors[0], cross_product)
- # Check if the numerator is nearly zero
- if np.abs(numerator) < 1e-16:
- raise ValueError("Vectors are nearly parallel")
-
- norms = [np.linalg.norm(v) for v in vectors]
-
- # Compute the denominator of the solid angle formula
- denominator = (
- norms[0] * norms[1] * norms[2] +
- np.dot(vectors[0], vectors[1]) * norms[2] +
- np.dot(vectors[0], vectors[2]) * norms[1] +
- np.dot(vectors[1], vectors[2]) * norms[0]
- )
-
- # Compute the solid angle in radians
- solid_angle = 2 * np.arctan2(numerator, denominator)
- # print("Solid angle:", solid_angle)
- return solid_angle
- except ZeroDivisionError:
- raise ValueError("Error computing solid angle: division by zero")
- except np.linalg.LinAlgError:
- raise ValueError("Error computing solid angle: linear algebra error")
-
-
- # BSE coefficient functions (Equations 11 and 10, Section 2.1.2)
- def eta_0(_A: float) -> float:
- """Equation 11 from the paper"""
- _eta_0: float = -0.0254 + 0.016 * _A - 1.86e-4 * _A ** 2 + 8.3e-7 * _A ** 3
- print('eta_0: ', _eta_0)
- return _eta_0
-
-
- # BSE Coefficient for single atomic targets, found to hold in the range of 10–100 keV
- def eta(theta: float, _A: float) -> float:
- """Equation 10 from the paper"""
- B: float = 0.89
- _eta: float = B * (eta_0(_A) / B) ** np.cos(theta)
- print('eta: ', _eta)
- return _eta
-
-
- '''
- # Compute the weighted average BSE coefficient for the alloy
- def compute_weighted_eta(theta: float, alloy: Alloy) -> float:
- """Equation 12 from the paper"""
- weighted_eta: float = 0
- for element, weight_fraction in zip(alloy.elements, alloy.composition):
- # print('element: ', element)
- # print('weight_fraction: ', weight_fraction)
- _A: float = element.atomic_number
- # print('A: ', _A)
- eta_i: float = eta(theta, _A)
- # print('eta_i: ', eta_i)
- weighted_eta += weight_fraction * eta_i
- # print('weighted_eta: ', weighted_eta)
- return weighted_eta
- '''
-
-
- # Correction function C(theta) (Equation 14, Section 2.1.3)
- def C(theta: float) -> float:
- print("theta1:", theta)
- """Equation 14 from the paper"""
- theta = np.rad2deg(theta) # Ensure angle is in degrees for C function
- _C: float = 1.01 - 4.06 * theta + 1.36e-5 * theta ** 2 - 1.71e-6 * theta ** 3
- print('C: ', _C)
- # Check if C is negative and raise an exception
- if _C < 0:
- raise ValueError(f"C function returned a negative value for theta: {theta}. "
- f"Please check the input and the C function implementation.")
- return _C
-
-
- # Scattering function g_BSE (Equation 13, Section 2.1.3)
- def g_BSE(_theta: float, _alpha: float, _beta: float) -> float:
- print("_theta, _alpha, _beta:", _theta, _alpha, _beta)
- """Equation 13 from the paper"""
- # k: float = np.rad2deg(_theta) / 13 # Correct implementation as per the paper
- # print("k: ", k)
- diffusive_part: float = (eta_0(A) / pi) * np.cos(_alpha)
- print(np.cos(_alpha))
- print("diffusive_part:", diffusive_part)
- # reflective_part: float = ((eta(_theta, A) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
- # if np.abs(reflective_part) < 1e-10:
- # reflective_part = 0
- # print("reflective_part:", reflective_part)
-
- # Compute a weighted average BSE coefficient for the alloy
- # reflective_part: float = ((compute_weighted_eta(_theta, Ti6Al4V) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
-
- # print('g_BSE: ', diffusive_part + reflective_part)
- return diffusive_part #+ reflective_part
-
-
- # Calculate the beam direction vector p_E
- def calculate_p_E(_beam_x, _beam_y):
- print(f"calculate_p_E({_beam_x}, {_beam_y}):")
- # Calculate the vector from the beam gun to the beam location
- beam_vector = np.array([_beam_x * dx, _beam_y * dy, -detector_height], dtype=np.float64)
- return normalize_vector(beam_vector)
-
- def calculate_p_E1(_beam_x, _beam_y, _height):
- """Calculate normalized beam direction vector from gun to surface point
-
- Args:
- _beam_x: Beam x-position relative to center (in cells)
- _beam_y: Beam y-position relative to center (in cells)
- _height: Height difference between gun and surface point (m)
- """
- print(f"calculate_p_E({_beam_x}, {_beam_y}):")
- print("_height:", _height)
-
- # Calculate vector components
- x_component = _beam_x * dx # Convert cells to meters
- y_component = _beam_y * dy # Convert cells to meters
- z_component = _height # Already in meters
-
- # Create beam vector (from gun to surface point)
- beam_vector = np.array([x_component, y_component, z_component], dtype=np.float64)
-
- # Verify downward direction
- if normalize_vector(beam_vector)[2] > 0:
- raise ValueError("Beam vector pointing upward")
-
- # Normalize and return
- return normalize_vector(beam_vector)
-
-
- # Calculate phi as the angle between p_E and the negative z-axis
- def calculate_phi(_p_E) -> float:
- z_axis = np.array([0.0, 0.0, -1.0], dtype=np.float64)
-
- # Check if _p_E is approximately equal to the negative z-axis
- '''if np.allclose(_p_E, z_axis, rtol=1e-8, atol=1e-8):
- _phi_rad = 0.0
- else:
- # Calculate the angle between _p_E and the z-axis using arccos
- _phi_rad = calculate_angle(_p_E, z_axis)'''
- _phi_rad = calculate_angle(_p_E, z_axis)
- print(f"calculate_phi({_p_E}, {_phi_rad}):")
- return _phi_rad
-
-
- def set_beam_movement(_x_center, _y_center, _x_offset, _y_offset, _movement_type):
- if _movement_type == 'X':
- if _x_offset >= 0:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset - 1
- else:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset + 1
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center + _y_offset + 1
- elif _movement_type == 'Y':
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center + _x_offset + 1
- if _y_offset >= 0:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset - 1
- else:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset + 1
- elif _movement_type == 'XY':
- if _x_offset >= 0:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset - 1
- else:
- _beam_x_start = _x_center + _x_offset
- _beam_x_end = _x_center - _x_offset + 1
- if _y_offset >= 0:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset - 1
- else:
- _beam_y_start = _y_center + _y_offset
- _beam_y_end = _y_center - _y_offset + 1
- else:
- raise ValueError("movement_type must be 'X' or 'Y' or 'XY'")
-
- return _beam_x_start, _beam_x_end, _beam_y_start, _beam_y_end
-
-
- # Initialize arrays to store electron counts
- electrons_per_detector_per_step = np.zeros((x_num, y_num, len(detector_positions)))
- total_electrons_per_step = np.zeros((x_num, y_num))
-
- # Additional arrays to capture detailed information
- electrons_per_detector_per_beam_profile = np.zeros(
- (beam_profile.shape[0], beam_profile.shape[1], len(detector_positions)))
- total_electrons_per_beam_profile = np.zeros((beam_profile.shape[0], beam_profile.shape[1]))
-
- # Calculate the center indices
- x_center = x_num // 2
- y_center = y_num // 2
- # print('x_center, y_center:', x_center, y_center)
-
- x_offset: int = 0
- y_offset: int = 0
- movement_type: str = 'Y'
-
- # Beam X-Movement
- # beam_x_start, beam_x_end = x_center - x_offset, x_center + x_offset + 1
- # beam_y_start, beam_y_end = y_center + y_offset, y_center + y_offset + 1
-
- # Beam Y-Movement
- # beam_x_start, beam_x_end = x_center + x_offset, x_center + x_offset + 1
- # beam_y_start, beam_y_end = y_center - y_offset, y_center + y_offset + 1
-
- # For X or Y Movement
- beam_x_start, beam_x_end, beam_y_start, beam_y_end = set_beam_movement(
- x_center, y_center, x_offset, y_offset, movement_type)
-
- # Determine the step size based on the direction of movement
- step_x = 1 if beam_x_start <= beam_x_end else -1
- step_y = 1 if beam_y_start <= beam_y_end else -1
-
- # print(beam_x_start, beam_x_end, beam_y_start, beam_y_end, step_x, step_y)
-
-
- z_heights = np.zeros(beam_size)
- fill_level_b = np.zeros(beam_size)
- fill_level_a = np.zeros(beam_size)
- normals = np.zeros((beam_size, 3))
- normals_b = np.zeros((beam_size, 3))
- normals_a = np.zeros((beam_size, 3))
- cell_counter = 0
-
-
- def normalize_vector(vector: np.ndarray) -> np.ndarray:
- """Normalize a 3D vector."""
- vector = vector.astype(np.float64) # Convert to float64
- norm = np.linalg.norm(vector)
- if norm < 1e-18: # Avoid division by zero
- return vector # Return original if nearly zero
- return vector / norm # Normalize
-
-
- def calculate_angle(_v1: np.ndarray, _v2: np.ndarray) -> float:
- """Calculate the angle between two normalized vectors."""
- _v1 = normalize_vector(_v1)
- _v2 = normalize_vector(_v2)
-
- dot_product = np.clip(np.dot(_v1, _v2), -1.0, 1.0)
-
- return np.arccos(dot_product) # Return angle in radians
-
-
- # Main simulation loop to iterate over the simulation domain
- for beam_x in range(beam_x_start, beam_x_end, step_x):
- for beam_y in range(beam_y_start, beam_y_end, step_y):
- beam_loc = (beam_x, beam_y)
- print(f"Beam center is at domain location: {beam_loc}")
-
- # Calculate the beam direction vector p_E
- # >> p_E = calculate_p_E(beam_x - x_center, beam_y - y_center)
- # >> print('p_E: ', p_E)
-
- # Calculate phi as the angle between p_E and the negative z-axis
- # phi_rad = np.arctan2(np.sqrt(dx**2 + dy**2), w_d)
- # >> phi_rad = calculate_phi(p_E)
- # print(f"phi_rad: {phi_rad}")
-
- # Reset counters for this beam location
- total_electrons: int = 0
- electrons_per_detector = np.zeros(len(detector_positions))
-
- # Iterate over the beam profile
- # for profile_x in [40]:
- for profile_x in range(beam_profile.shape[0]):
- # for profile_y in [79]:
- for profile_y in range(beam_profile.shape[1]):
- print(f"Beam profile shape ({profile_x}, {profile_y})")
- domain_x = beam_loc[0] - beam_profile.shape[0] // 2 + profile_x
- domain_y = beam_loc[1] - beam_profile.shape[1] // 2 + profile_y
- print(f"Domain ({domain_x}, {domain_y})")
-
- # p_E = calculate_p_E(domain_x - x_center, domain_y - y_center)
- # print('p_E: ', p_E)
- # phi_rad = calculate_phi(p_E)
-
- if (0 <= domain_x < x_num) and (0 <= domain_y < y_num):
- if beam_profile[profile_x, profile_y] > threshold:
- print('finding interface_z_indices')
- # print('interface_z_indices:', interface_z)
- z_interp, z_height, fill_level_below, fill_level_above, normal_interp, normals_below, normals_above \
- = interpolate_fill_level_and_normal_new(fill_levels_smoothed, n_S, domain_y, domain_x)
- if z_interp is not None:
- z_heights[cell_counter] = z_height
- fill_level_b[cell_counter] = fill_level_below
- fill_level_a[cell_counter] = fill_level_above
- normals[cell_counter] = normal_interp
- normals_b[cell_counter] = normals_below
- normals_a[cell_counter] = normals_above
- cell_counter += 1
- n_S_local = normal_interp
- print('n_S_local:', n_S_local)
- if not np.allclose(np.linalg.norm(n_S_local), 1.0, atol=1e-5):
- raise ValueError(f"n_S_local is not normalized, {np.linalg.norm(n_S_local)}")
-
- p_E = calculate_p_E1(domain_x - x_center, domain_y - y_center, -(detector_height-z_height))
- print('p_E: ', p_E)
-
- theta_local = calculate_angle(-p_E, n_S_local)
- print('theta_local:', theta_local)
- # if theta_local < 0.0 or theta_local > 0.011:
- # raise ValueError(f"theta_local value should be 0.0, found {theta_local}")
- # theta_local: 0.016
-
- scattering_point = np.array([
- (domain_x - beam_x) * dx,
- (domain_y - beam_y) * dy,
- z_height
- ], dtype=np.float64)
- print(f'scattering_point ({domain_x}, {domain_y}, {z_interp}):', scattering_point)
-
- for det_idx, det_pos in enumerate(detector_positions):
- # print(f"Detector {det_idx + 1} position: {det_pos}")
- det_idx: int
- # det_pos: np.ndarray[float]
- print(f'det_pos[{det_idx + 1}]: {det_pos}')
-
- # Calculate the vector from the scattering point to the detector
- det_vec = det_pos - scattering_point
- det_vec = normalize_vector(det_vec)
- print(f'det_vec[{det_idx + 1}]: {det_vec}')
- if not np.allclose(np.linalg.norm(det_vec), 1.0):
- raise ValueError("det_vec is not normalized.")
-
- # alpha = np.arccos(np.clip(np.dot(det_vec, n_S_local), -1.0, 1.0))
- alpha = calculate_angle(det_vec, n_S_local)
- print(f'alpha: {alpha}')
-
- if alpha > np.pi / 2:
- continue
-
- p_E_reflected = p_E - 2.0 * np.dot(p_E, n_S_local) * n_S_local
- # p_E_reflected = np.where(np.abs(p_E_reflected) < 1e-8, 0.0, p_E_reflected)
- p_E_reflected = normalize_vector(p_E_reflected)
- print(f'p_E_reflected: {p_E_reflected}')
-
- # beta = np.arccos(np.clip(np.dot(p_E_reflected, det_vec), -1.0, 1.0))
- beta = calculate_angle(p_E_reflected, det_vec)
- print(f'beta: {beta}')
-
- # Apply energy threshold
- # electron_energy = N_0 * e * g_bse_value
- # print('electron_energy: ', electron_energy)
-
- # Calculate the initial energy of the primary electrons (E_0)
- V_acc = 60e3 # Accelerating voltage in Volts
- E_0 = V_acc * e # Initial energy of primary electrons in joules
-
- # Calculate the energy of the backscatter electrons
- E_BSE = E_0 * eta(theta_local, A)
- # Compute a weighted average BSE coefficient for the alloy
- # E_BSE = E_0 * compute_weighted_eta(theta_local, Ti6Al4V)
- # print('E_BSE: ', E_BSE)
-
- # Bifurcation logic for PE and SE
- # if E_BSE > threshold_energy:
- # Backscatter electrons (BSE)
- # Compute scattering function g_bse_value
- g_bse_value = g_BSE(theta_local, alpha, beta)
- print(f'g_bse_value: {g_bse_value}')
-
- # Ensure g_bse_value is non-negative
- if g_bse_value <= 0:
- raise ValueError(f"g_bse_value is {g_bse_value} (zero negative), check calculations.")
-
- # Compute solid angle (Equation A1)
- v1, v2, v3 = get_detector_vertices(det_pos)
- normal = np.cross(v2 - v1, v3 - v1)
- normal = normalize_vector(normal)
-
- dOmega = compute_solid_angle(v1 - scattering_point,
- v2 - scattering_point,
- v3 - scattering_point)
- print(f"Detector {det_idx + 1} solid angle: {dOmega}")
-
- # Compute number of electrons collected by the detector (Equation A2)
- # The beam ray stays on each grid point ij for the dwell time t_d and distributes N_0 PE into it
- N_0 = beam_current * dwell_time / e # Equation 3, Section 2.1.1
- print("N_0:", N_0)
-
- # There are N1 electrons in the chamber after the 1st scattering event
- N_1 = N_0 * eta(theta_local, A) # Equation 4, Section 2.1.1
- # Compute a weighted average BSE coefficient for the alloy
- # N_1 = N_0 * compute_weighted_eta(theta_local, Ti6Al4V) # Equation 4, Section 2.1.1
- print("N_1:", N_1)
-
- # The number of electrons N_1_d collected by a general detector d
- N_1_d = N_1 * np.sum(g_bse_value * dOmega) # Equation 5, Section 2.1.1
-
- # Accumulate electrons per detector for this beam profile point
- electrons_per_detector_per_beam_profile[profile_x, profile_y, det_idx] += N_1_d
-
- # Accumulate electrons per detector for this domain point
- electrons_per_detector_per_step[beam_x, beam_y, det_idx] += N_1_d
-
- # Add to total electrons for this beam profile point
- total_electrons_per_beam_profile[profile_x, profile_y] += N_1_d
-
- # Add to total electrons for this domain point
- total_electrons_per_step[beam_x, beam_y] += N_1_d
-
- # Accumulate electrons per detector
- electrons_per_detector[det_idx] += N_1_d
-
- # Add to total electrons
- total_electrons += N_1_d
-
- # print(f"After addition: electrons_per_detector[{det_idx + 1}]: {electrons_per_detector[det_idx]}")
- else:
- print(f"No interface found at ({domain_x}, {domain_y})")
-
- # Print detailed information for this beam profile point
- # print(f"Electrons per detector at this beam profile point ({profile_x}, {profile_y}): {electrons_per_detector_per_beam_profile[profile_x, profile_y]}")
- # print(f"Total electrons at this beam profile point ({profile_x}, {profile_y}): {total_electrons_per_beam_profile[profile_x, profile_y]}")
-
- # Print detailed information for this domain point after iterating over the beam profile
- # print(f"Electrons per detector at this domain point ({beam_x}, {beam_y}): {electrons_per_detector_per_step[beam_x, beam_y]}")
- # print(f"Total electrons at this domain point ({beam_x}, {beam_y}): {total_electrons_per_step[beam_x, beam_y]}")
-
- # print("Shape of electrons_per_detector_per_step:", electrons_per_detector_per_step.shape)
- # print("Number of lines:", len(lines))
-
- # Compute final results
- total_electrons_all_steps = np.sum(total_electrons_per_step)
- total_electrons_per_detector = np.sum(electrons_per_detector_per_step, axis=(0, 1))
-
- # Print final results
- # print(f"Total electrons hitting all detectors: {total_electrons_all_steps}")
- # for i, electrons in enumerate(total_electrons_per_detector):
- # print(f"Total electrons hitting detector {i + 1}: {electrons}")
-
- # for beam_x in range(beam_x_start, beam_x_end, step_x):
- # for beam_y in range(beam_y_start, beam_y_end, step_y):
- # print(f"Electrons per detector at this domain point ({beam_x}, {beam_y}): {electrons_per_detector_per_step[beam_x, beam_y]}")
- # print('end simulation')
- # print(electrons_per_detector_per_step)
-
- # print("total_electrons:", total_electrons, "at timestep", current_timestep)
- # print("electrons_per_det:", electrons_per_det, "at timestep", current_timestep)
- # print("electrons_per_beam_loc:", electrons_per_beam_loc, "at timestep", current_timestep)
- # print("electrons_per_beam_loc_per_det:", electrons_per_beam_loc_per_det, "at timestep", current_timestep)
- return electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, normals, normals_b, normals_a
-
-def read_fill_levels(filepath):
- """Read only the FillingFrac grid from VDB file"""
- try:
- # Try to read just the FillingFrac grid
- grid = vdb.read(filepath, "FillingFrac")
-
- # Get dimensions
- dims = grid.evalActiveVoxelDim()
- print(f"Grid dimensions: {dims}")
-
- # Create array
- fill_array = np.zeros(dims, dtype=np.float64)
-
- # Copy data
- grid.copyToArray(fill_array)
-
- # Transpose to ray tracer's expected order (z, y, x)
- fill_array = np.transpose(fill_array, (2, 1, 0))
-
- # Create standardized array with 256 z cells
- standardized_array = np.zeros((256, dims[1], dims[0]), dtype=np.float64)
-
- # Copy data from original array to standardized array
- standardized_array[:fill_array.shape[0], :, :] = fill_array
-
- # Apply Gaussian smoothing
- sigma = 1.0
- standardized_array_smoothed = ndimage.gaussian_filter(standardized_array, sigma=sigma, mode='nearest', radius=1)
-
- # Print values
- x = 157
- y = 196
- print(f"\nFill levels at ({x}, {y}):")
- print("z_index | fill_level")
- print("-" * 20)
-
- # Get fill levels
- fill_levels = standardized_array_smoothed[:, y, x]
- valid_levels = fill_levels[(fill_levels > 0) & (fill_levels < 1)]
-
- if len(valid_levels) > 0:
- for z, val in enumerate(fill_levels):
- if 0 <= val <= 1:
- print(f"{z:7d} | {val:.6f}")
- else:
- print("No fill level found between 0 and 1")
-
- return standardized_array_smoothed
-
- except Exception as e:
- print(f"Error reading {filepath}")
- print(f"Error details: {str(e)}")
- try:
- metadata = vdb.readMetadata(filepath)
- print("\nFile metadata:")
- for key, value in metadata.items():
- print(f"{key}: {value}")
- except Exception as meta_e:
- print(f"Error reading metadata: {str(meta_e)}")
- return None
-
-def process_single_timestep(filepath, timestep):
- """Process a single VDB file and run ray tracer"""
- try:
- # Read the FillingFrac grid
- grid = vdb.read(filepath, "FillingFrac")
-
- # Get dimensions
- dims = grid.evalActiveVoxelDim()
- x_num, y_num = dims[0], dims[1]
- z_num = 256 # dims[2] # Fixed z dimension
- print(f"Domain dimensions: {x_num} x {y_num} x {z_num}")
-
- # Read fill levels
- fill_levels = read_fill_levels(filepath)
- if fill_levels is None:
- raise ValueError("Could not read fill levels from VDB file")
-
- if fill_levels is not None:
- print("\nSuccessfully read fill levels")
- print(f"Array shape: {np.shape(fill_levels)}")
-
- try:
- # Convert to numpy array if not already
- fill_array = np.asarray(fill_levels)
-
- # Get non-zero elements
- non_zero = fill_array[fill_array > 0]
-
- # Get interface cells (0 < fill_level < 1)
- interface_cells = fill_array[(fill_array > 0) & (fill_array < 1)]
-
- # Print statistics
- if len(non_zero) > 0:
- print("\nNon-zero cells statistics:")
- print(f"Number of non-zero cells: {len(non_zero)}")
- print(f"Min non-zero value: {np.min(non_zero):.6f}")
- print(f"Max value: {np.max(fill_array):.6f}")
-
- if len(interface_cells) > 0:
- print("\nInterface cells statistics:")
- interface_percentage = (len(interface_cells) * 100.0) / len(non_zero)
- print(f"Number of interface cells: {len(interface_cells)} ({interface_percentage:.3f}%)")
- print(f"Min interface value: {np.min(interface_cells):.6f}")
- print(f"Max interface value: {np.max(interface_cells):.6f}")
- print(f"Mean interface value: {np.mean(interface_cells):.6f}")
- else:
- print("\nNo interface cells found (no fill levels between 0 and 1)")
-
- except Exception as e:
- print(f"Error processing array: {str(e)}")
-
- # Run ray tracer
- print(f"\nProcessing timestep {timestep}")
- electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, \
- normals, normals_b, normals_a = runRayTracer(
- x_num,
- y_num,
- z_num,
- # fill_levels.flatten(),
- fill_levels,
- timestep
- )
-
- return electrons_per_detector_per_step, z_heights, fill_level_b, fill_level_a, \
- normals, normals_b, normals_a
-
- except Exception as e:
- print(f"Error processing timestep {timestep}: {str(e)}")
- raise
-
-def main():
- """Main function to process VDB files and run ray tracer"""
- # Directory containing VDB files
- directory = "/local/ca00xebo/softwares/KiSSAM/Markl-plate"
-
- # Get sorted list of VDB files
- vdb_files = []
- for f in os.listdir(directory):
- if f.startswith('geometry-00000') and f.endswith('.vdb'):
- try:
- # Extract timestep number from filename
- timestep = int(f.split('-')[1].split('.')[0])
- vdb_files.append((timestep, f))
- except (IndexError, ValueError) as e:
- print(f"Skipping file {f}: {str(e)}")
-
- # Sort by timestep
- vdb_files.sort() # sort based on timestep number
-
- total_timesteps = len(vdb_files)
- print(f"Found {total_timesteps} VDB files to process")
- electrons_across_timesteps = np.zeros((total_timesteps, 4))
-
- # Process each timestep
- _results = []
- for timestep, vdb_file in enumerate(vdb_files):
- filepath = os.path.join(directory, vdb_file[1]) # vdb_file[1] contains the filename
- try:
- print(f"\nProcessing file {vdb_file[1]} (timestep {timestep})")
-
- # Process this timestep
- _result = process_single_timestep(filepath, timestep)
-
- # Debug electron counts
- print("\nElectron count debug:")
- print(f"Shape of electrons_per_detector_per_step: {_result[0].shape}")
- print(f"Min value: {np.min(_result[0])}")
- print(f"Max value: {np.max(_result[0])}")
- print(f"Mean value: {np.mean(_result[0])}")
-
- # Store electron counts with debugging
- total_electrons_for_detectors = np.sum(_result[0], axis=(0, 1))
- print(f"Total electrons before assignment: {total_electrons_for_detectors}")
-
- electrons_across_timesteps[timestep] = total_electrons_for_detectors
- print(f"Electrons for timestep {timestep}:")
- print(f"Raw values: {electrons_across_timesteps[timestep]}")
- print(f"Sum: {np.sum(electrons_across_timesteps[timestep])}")
- print("electrons_across_timesteps:", electrons_across_timesteps)
-
- _results.append(_result)
- print(f"Completed timestep {timestep}")
-
- except Exception as e:
- print(f"Error processing timestep {timestep}: {str(e)}")
- raise
-
- # Final debug output
- print("\nFinal electrons_across_timesteps:")
- print(f"Shape: {electrons_across_timesteps.shape}")
- print(f"Min value: {np.min(electrons_across_timesteps)}")
- print(f"Max value: {np.max(electrons_across_timesteps)}")
- print(f"Mean value: {np.mean(electrons_across_timesteps)}")
-
- # Plot results
- plot_detector_results1(electrons_across_timesteps)
-
- return _results
-
-def plot_detector_results(electrons_per_detector_time):
- """Plot electron counts for each detector across timesteps"""
- plt.figure(figsize=(12, 8))
- timesteps = range(len(electrons_per_detector_time))
-
- # Plot for each detector
- for detector in range(4):
- plt.plot(timesteps, electrons_per_detector_time[:, detector],
- label=f'Detector {detector+1}',
- marker='o')
-
- plt.title("Electrons Captured by Detectors Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Number of Electrons")
- plt.legend()
- plt.grid(True)
-
- # Add value annotations
- for detector in range(4):
- for t, value in zip(timesteps, electrons_per_detector_time[:, detector]):
- plt.annotate(f'{value:.2e}',
- (t, value),
- textcoords="offset points",
- xytext=(0,10),
- ha='center',
- rotation=45)
-
- plt.tight_layout()
- plt.savefig("detector_electrons_over_time(VDB-multi).png", dpi=300, bbox_inches='tight')
- plt.close()
-
- # Create additional visualization - normalized values
- plt.figure(figsize=(12, 8))
- normalized_electrons = electrons_per_detector_time / np.max(electrons_per_detector_time)
-
- for detector in range(4):
- plt.plot(timesteps, normalized_electrons[:, detector],
- label=f'Detector {detector+1}',
- marker='o')
-
- plt.title("Normalized Electrons Captured by Detectors Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Normalized Electron Count")
- plt.legend()
- plt.grid(True)
- plt.tight_layout()
- plt.savefig("detector_electrons_over_time_normalized(VDB-multi).png", dpi=300, bbox_inches='tight')
- plt.close()
-
-def plot_detector_results1(electrons_per_detector_time):
- """
- Plot electron counts for detector pairs across timesteps
-
- Args:
- electrons_per_detector_time: Array of shape (num_timesteps, 4) containing
- electron counts for each detector at each timestep
- :param electrons_per_detector_time:
- """
- # Create figure for detectors 1 and 2
- plt.figure(figsize=(12, 8))
- timesteps = range(len(electrons_per_detector_time))
-
- # Plot detectors 1 and 2
- plt.plot(timesteps, electrons_per_detector_time[:, 0],
- label='Detector 1', marker='o', linewidth=1, markevery=20)
- plt.plot(timesteps, electrons_per_detector_time[:, 1],
- label='Detector 2', marker='s', linewidth=1, markevery=20)
-
- plt.title("Electrons Captured by Detectors 1 and 2 Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Number of Electrons")
- plt.legend()
- plt.grid(True)
- plt.tight_layout()
- plt.savefig("detectors_1_2_over_time(VDB-multi).png", dpi=300, bbox_inches='tight')
- plt.close()
-
- # Create figure for detectors 3 and 4
- plt.figure(figsize=(12, 8))
-
- # Plot detectors 3 and 4
- plt.plot(timesteps, electrons_per_detector_time[:, 2],
- label='Detector 3', marker='o', linewidth=1, markevery=20)
- plt.plot(timesteps, electrons_per_detector_time[:, 3],
- label='Detector 4', marker='s', linewidth=1, markevery=20)
-
- plt.title("Electrons Captured by Detectors 3 and 4 Over Time")
- plt.xlabel("Timestep")
- plt.ylabel("Number of Electrons")
- plt.legend()
- plt.grid(True)
- plt.tight_layout()
- plt.savefig("detectors_3_4_over_time(VDB-multi).png", dpi=300, bbox_inches='tight')
- plt.close()
-
-if __name__ == "__main__":
- results = main()
\ No newline at end of file
diff --git a/apps/showcases/FreeSurface/raytracerVDB.py b/apps/showcases/FreeSurface/raytracerVDB.py
index 6ff1874304a89a2d6b7eeab0239cf072fdf6a138..e10db3bff59e730fb7f9d69ed10d0b888bfd8dcf 100644
--- a/apps/showcases/FreeSurface/raytracerVDB.py
+++ b/apps/showcases/FreeSurface/raytracerVDB.py
@@ -1555,6 +1555,9 @@ def main():
# Sort by timestep
vdb_files.sort() # sort based on timestep number
+ # Slice the first 120 files
+ # vdb_files = vdb_files[:120]
+
total_timesteps = len(vdb_files)
print(f"Found {total_timesteps} VDB files to process")
electrons_across_timesteps = np.zeros((total_timesteps, 4))
diff --git a/apps/showcases/FreeSurface/testRayTracerFunctions.py b/apps/showcases/FreeSurface/testRayTracerFunctions.py
deleted file mode 100644
index 58afe1f4d71fd8ad90a846177180800f042d4019..0000000000000000000000000000000000000000
--- a/apps/showcases/FreeSurface/testRayTracerFunctions.py
+++ /dev/null
@@ -1,321 +0,0 @@
-import pyopenvdb as vdb
-import os
-from math import floor
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib.colors as colors
-from matplotlib.animation import FuncAnimation
-from typing import Tuple, List
-import scipy.ndimage as ndimage
-
-def runRayTracer(x_num, y_num, z_num, fill_levels, current_timestep):
-
- print("current_timestep", current_timestep)
- # fill_levels = np.clip(fill_levels, 0, 1)
- fill_levels = fill_levels.astype(np.float64)
-
- # Define the tolerance
- tolerance = 1e-8
-
- # Check if any value in fill_levels is outside the [0, 1] range
- out_of_range = np.logical_or(fill_levels < -tolerance, fill_levels > 1 + tolerance)
- if np.any(out_of_range):
- out_of_range_values = fill_levels[out_of_range]
- out_of_range_indices = np.where(out_of_range)
-
- for idx, value in zip(zip(*out_of_range_indices), out_of_range_values):
- print(f"Index: {idx}, Value: {value}")
-
- # raise ValueError(f"Fill level values must be between 0 and 1 inclusive (with tolerance of {tolerance}).")
-
- # Check for fill levels within the range [0, 1]
- in_range = np.logical_and(fill_levels > 0, fill_levels < 1)
- in_range_values = fill_levels[in_range]
- in_range_indices = np.where(in_range)
-
- # print("Fill level values within range (0 to 1):")
- # for idx, value in zip(zip(*in_range_indices), in_range_values):
- # if idx == (np.int64(309),):
- # print(f"Index: {idx}, Value: {value}")
-
- # Constants
- pi: float = np.pi # Mathematical constant π
- e: float = 1.60217662e-19 # Elementary charge in Coulombs
- beam_current: float = 3e-3 # Beam current in Amperes
- dwell_time: float = 0.32e-6 # Dwell time in seconds
- A: float = 21.5 # Approximation for Ti6Al4V
- # Ti6Al4V = Ti6Al4V.create_Ti6Al4V(T)
- # A: float = Ti6Al4V.atomic_number
- # print('Ti6Al4V.elements: ', Ti6Al4V.elements)
- # print('Ti6Al4V.composition: ', Ti6Al4V.composition)
- # print('A: ', A)
- dx: float = 3e-6 # Cell size in meters
- dy: float = dx # Cell size in meters
- dz: float = dx # Cell size in meters
- threshold: float = 0 # Beam profile threshold
- # The voltage at detector d can be computed using the transimpedance T (in units of V /A ) of a virtual amplifier in the context of the model
- transimpedance: float = 1e5 # Transimpedance in Volts per Ampere
- eV_threshold = 50 # Threshold energy in eV
- threshold_energy = eV_threshold * e # Threshold energy in joules
-
- # Domain setup (Section 2.1.1)
- # Define the 3D domain for the simulation
- x_min, x_max = -0.5 * x_num * dx, 0.5 * x_num * dx # -10 mm to 10 mm
- # print('x_min, x_max:', x_min, x_max)
- y_min, y_max = -0.5 * y_num * dy, 0.5 * y_num * dy # -10 mm to 10 mm
- # print('y_min, y_max:', y_min, y_max)
- z_min, z_max = 0, z_num * dz
- # print('z_min, z_max:', z_min, z_max)
-
- fill_levels_reshaped = np.reshape(fill_levels, (z_num, y_num, x_num))
- # fill_levels_reshaped = fill_levels
- # print('fill_levels_reshaped:', fill_levels_reshaped)
- out_of_range_fill_levels = np.logical_or(fill_levels_reshaped < -tolerance, fill_levels_reshaped > 1 + tolerance)
- if np.any(out_of_range_fill_levels):
- out_of_range_fill_level_values = fill_levels_reshaped[out_of_range_fill_levels]
- out_of_range_fill_level_indices = np.where(out_of_range_fill_levels)
-
- for idx, value in zip(zip(*out_of_range_fill_level_indices), out_of_range_fill_level_values):
- print(f"Index: {idx}, Value: {value}")
- raise ValueError(f"Fill level values must be between 0 and 1 inclusive (with tolerance of {tolerance}).")
-
- sigma = 0.0
- # Different sigma for each dimension
- # sigma = [0.0, 1.0, 0.5] # [z, y, x]
- fill_levels_smoothed = ndimage.gaussian_filter(fill_levels_reshaped, sigma=sigma, mode='nearest', radius=1)
-
- fill_levels_smoothed = np.clip(
- fill_levels_smoothed,
- a_min=np.float64(0.0),
- a_max=np.float64(1.0)
- )
-
- # Check if values are still in bounds
- if not (np.all(fill_levels_smoothed >= 0) and np.all(fill_levels_smoothed <= 1)):
- raise ValueError("Smoothed fill levels out of bounds!")
-
- def normalize_vector(vector: np.ndarray) -> np.ndarray:
- """Normalize a 3D vector."""
- vector = vector.astype(np.float64) # Convert to float64
- norm = np.linalg.norm(vector)
- if norm < 1e-18: # Avoid division by zero
- return vector # Return original if nearly zero
- return vector / norm # Normalize
-
- def interpolate_normals_new(_normal_below: np.ndarray, _normal_above: np.ndarray, t: float) -> np.ndarray:
- """Interpolate between two normals using SLERP"""
- if not (0 <= t <= 1):
- raise ValueError(f"Interpolation parameter t must be in [0,1], got {t}")
-
- # Ensure normalized vectors
- _normal_below = normalize_vector(_normal_below.astype(np.float64))
- _normal_above = normalize_vector(_normal_above.astype(np.float64))
-
- # Compute dot product
- dot_product = np.clip(np.dot(_normal_below, _normal_above), -1.0, 1.0)
-
- # Thresholds for special cases
- PARALLEL_THRESHOLD = 0.9999
- ANTIPARALLEL_THRESHOLD = -0.9999
-
- # Case 1: Nearly parallel
- if dot_product > PARALLEL_THRESHOLD:
- return normalize_vector((1.0 - t) * _normal_below + t * _normal_above)
-
- # Case 2: Nearly antiparallel
- if dot_product < ANTIPARALLEL_THRESHOLD:
- perp = np.array([1.0, 0.0, 0.0])
- if abs(np.dot(_normal_below, perp)) > 0.9:
- perp = np.array([0.0, 1.0, 0.0])
- perp = normalize_vector(np.cross(_normal_below, perp))
- angle = np.pi * t
- return normalize_vector(_normal_below * np.cos(angle) + perp * np.sin(angle))
-
- # Case 3: Standard SLERP
- theta = np.arccos(dot_product)
- sin_theta = np.sin(theta)
-
- if abs(sin_theta) < 1e-10:
- return normalize_vector(_normal_below)
-
- scale0 = np.sin((1.0 - t) * theta) / sin_theta
- scale1 = np.sin(t * theta) / sin_theta
- return normalize_vector(scale0 * _normal_below + scale1 * _normal_above)
-
-
- def interpolate_fill_level_and_normal_new(_fill_level, surface_normals, _y, _x):
- """Find interface points and interpolate normals"""
- if surface_normals.ndim != 4:
- raise ValueError(f"Incorrect shape for surface normals: {surface_normals.ndim}")
-
- _column = _fill_level[:, _y, _x]
-
- # Constants
- INTERFACE_THRESHOLD = 0.5
- GRADIENT_THRESHOLD = 0.1
- TOLERANCE = 1e-6
-
- # Scan from top to bottom
- for z in range(len(_column)-2, -1, -1):
- current_fill = _column[z]
- next_fill = _column[z + 1]
-
- gradient = abs(next_fill - current_fill)
- if gradient > GRADIENT_THRESHOLD:
- # Use tolerance for interface detection
- crosses_interface = (
- (current_fill >= INTERFACE_THRESHOLD - TOLERANCE and
- next_fill <= INTERFACE_THRESHOLD + TOLERANCE)
- )
-
- if crosses_interface:
- print(f"Found interface at z={z}")
- print(f"fill_level_below: {current_fill}")
- print(f"fill_level_above: {next_fill}")
-
- # Linear interpolation
- denominator = next_fill - current_fill
- if abs(denominator) < TOLERANCE:
- _z_interp = z + 0.5
- t_interface = 0.5
- else:
- # Handle numerical precision near 0.5
- if abs(current_fill - INTERFACE_THRESHOLD) < TOLERANCE:
- t_interface = 0.0
- elif abs(next_fill - INTERFACE_THRESHOLD) < TOLERANCE:
- t_interface = 1.0
- else:
- t_interface = (INTERFACE_THRESHOLD - current_fill) / denominator
-
- # Ensure t is in valid range
- if not (0 <= t_interface <= 1):
- # Try adjusting the threshold slightly
- adjusted_threshold = INTERFACE_THRESHOLD
- if t_interface < 0:
- adjusted_threshold -= TOLERANCE
- else:
- adjusted_threshold += TOLERANCE
-
- t_interface = (adjusted_threshold - current_fill) / denominator
-
- if not (0 <= t_interface <= 1):
- error_msg = (
- f"\nInterpolation parameter t_interface={t_interface:.6f} is out of bounds [0,1]\n"
- f"Debug information:\n"
- f" - z-index: {z}\n"
- f" - current_fill: {current_fill:.6f}\n"
- f" - next_fill: {next_fill:.6f}\n"
- f" - denominator: {denominator:.6f}\n"
- f" - INTERFACE_THRESHOLD: {INTERFACE_THRESHOLD}\n"
- f" - Position (y,x): ({_y}, {_x})\n"
- f" - Fill level values around interface:\n"
- )
- for i in range(max(0, z-2), min(len(_column), z+3)):
- error_msg += f" {i:3d} | {_column[i]:.6f}\n"
- raise ValueError(error_msg)
-
- _z_interp = z + t_interface
-
- _z_height = (_z_interp + 0.5) * dz
-
- # Get and interpolate normals
- _normal_below = surface_normals[z, _y, _x].astype(np.float64)
- _normal_above = surface_normals[z + 1, _y, _x].astype(np.float64)
-
- _normal_interp = interpolate_normals_new(_normal_below, _normal_above, t_interface)
-
- return _z_interp, _z_height, next_fill, current_fill, \
- _normal_interp, _normal_below, _normal_above
-
- return None, None, None, None, None, None, None
-
-
- def compute_surface_normals_sobel_zeros(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with default upward normal
- _n_S_normalized = np.zeros_like(_n_S)
- _n_S_normalized[..., 2] = 1.0
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
-
- def compute_surface_normals_sobel_nans(_fill_levels: np.ndarray) -> np.ndarray:
- """Compute surface normals using Sobel filter with NaN initialization"""
- _fill_levels = _fill_levels.astype(np.float64)
-
- if _fill_levels.ndim != 3:
- raise ValueError("Input must be a 3D array")
-
- TOLERANCE = {
- 'interface': 1e-10,
- 'norm': 1e-10
- }
-
- # Compute gradients without initial inversion
- _grad_x = ndimage.sobel(_fill_levels, axis=2, output=np.float64, mode='nearest')
- _grad_y = ndimage.sobel(_fill_levels, axis=1, output=np.float64, mode='nearest')
- _grad_z = ndimage.sobel(_fill_levels, axis=0, output=np.float64, mode='nearest')
-
- # Stack gradients
- _n_S = np.stack([_grad_x, _grad_y, _grad_z], axis=-1)
-
- # Create interface mask
- interface_mask = np.logical_and(
- _fill_levels > TOLERANCE['interface'],
- _fill_levels < 1-TOLERANCE['interface']
- )[..., np.newaxis]
-
- # Ensure upward-pointing normals (from liquid to gas) after gradient calculation
- mask = _n_S[..., 2] < 0
- _n_S[mask] = -_n_S[mask]
-
- # Normalize
- norm = np.sqrt(np.sum(_n_S * _n_S, axis=-1, keepdims=True))
- norm = np.maximum(norm, TOLERANCE['norm'])
-
- # Initialize with NaN
- _n_S_normalized = np.full_like(_n_S, np.nan)
-
- # Apply normalization only at interface
- _n_S_normalized = np.where(interface_mask, _n_S / norm, _n_S_normalized)
-
- return _n_S_normalized
-
-
- # Compute surface normals (Section 2.1.1)
- n_S = compute_surface_normals_sobel_nans(fill_levels_smoothed)