Skip to content

generate.disk

Elastic disk solution and photoelastic simulation.

This module provides functions for simulating the photoelastic response of an elastic disk under compression, including analytical stress solutions and Mueller matrix-based polarimetry simulation.

disk

Functions

diametrical_stress_cartesian(X, Y, P, R)

Exact Brazil test solution from ISRM standards and Jaeger & Cook P: total load (force per unit thickness) R: disk radius

Key validation: At center (0,0): - sigma_x = 2P/(piR) (tensile) - sigma_y = -6P/(piR) (compressive) - tau_xy = 0

Source code in photoelastimetry/generate/disk.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def diametrical_stress_cartesian(X, Y, P, R):
    """
    Exact Brazil test solution from ISRM standards and Jaeger & Cook
    P: total load (force per unit thickness)
    R: disk radius

    Key validation: At center (0,0):
    - sigma_x = 2P/(pi*R) (tensile)
    - sigma_y = -6P/(pi*R) (compressive)
    - tau_xy = 0
    """

    X_safe = X.copy()
    Y_safe = Y.copy()

    # Small offset to avoid singularities at origin
    origin_mask = (X**2 + Y**2) < (0.001 * R) ** 2
    X_safe = np.where(origin_mask, 0.001 * R, X_safe)
    Y_safe = np.where(origin_mask, 0.001 * R, Y_safe)

    # Distance from load points
    r1 = np.sqrt(X_safe**2 + (Y_safe - R) ** 2)  # from (0, R)
    r2 = np.sqrt(X_safe**2 + (Y_safe + R) ** 2)  # from (0, -R)

    # Angles from load points
    theta1 = np.arctan2(X_safe, Y_safe - R)
    theta2 = np.arctan2(X_safe, Y_safe + R)

    sigma_xx = (
        -(2 * P / np.pi)
        * (np.cos(theta1) ** 2 * (Y_safe - R) / (r1**2) - np.cos(theta2) ** 2 * (Y_safe + R) / (r2**2))
        / R
    )

    sigma_yy = (
        -(2 * P / np.pi)
        * (np.sin(theta1) ** 2 * (Y_safe - R) / (r1**2) - np.sin(theta2) ** 2 * (Y_safe + R) / (r2**2))
        / R
    )

    tau_xy = (
        -(2 * P / np.pi)
        * (
            np.sin(theta1) * np.cos(theta1) * (Y_safe - R) / (r1**2)
            - np.sin(theta2) * np.cos(theta2) * (Y_safe + R) / (r2**2)
        )
        / R
    )

    return sigma_xx, sigma_yy, tau_xy

generate_synthetic_brazil_test(X, Y, P, R, S_i_hat, mask, wavelengths_nm, thickness, C, polarisation_efficiency)

Generate synthetic Brazil test data for validation This function creates a synthetic dataset based on the analytical solution and saves it in a format suitable for testing.

Source code in photoelastimetry/generate/disk.py 
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
def generate_synthetic_brazil_test(
    X, Y, P, R, S_i_hat, mask, wavelengths_nm, thickness, C, polarisation_efficiency
):
    """
    Generate synthetic Brazil test data for validation
    This function creates a synthetic dataset based on the analytical solution
    and saves it in a format suitable for testing.
    """

    # Get stress components directly
    sigma_xx, sigma_yy, tau_xy = diametrical_stress_cartesian(X, Y, P, R)

    # Mask outside the disk
    sigma_xx[~mask] = np.nan
    sigma_yy[~mask] = np.nan
    tau_xy[~mask] = np.nan

    # Principal stress difference and angle
    sigma_avg = 0.5 * (sigma_xx + sigma_yy)
    R_mohr = np.sqrt(((sigma_xx - sigma_yy) / 2) ** 2 + tau_xy**2)
    sigma1 = sigma_avg + R_mohr
    sigma2 = sigma_avg - R_mohr
    principal_diff = sigma1 - sigma2
    theta_p = 0.5 * np.arctan2(2 * tau_xy, sigma_xx - sigma_yy)

    # Mask again
    principal_diff[~mask] = np.nan
    theta_p[~mask] = np.nan

    height, width = sigma_xx.shape

    synthetic_images = np.empty((height, width, 3, 4))  # RGB, 4 polarizer angles

    # Use incoming light fully S1 polarized (standard setup)
    # S_i_hat = np.array([0.0, 0.0, 1.0])
    nu = 1.0  # Solid sample

    for i, lambda_light in tqdm(enumerate(wavelengths_nm)):
        # Generate four-step polarimetry images using Mueller matrix approach
        I0_pol, I45_pol, I90_pol, I135_pol = simulate_four_step_polarimetry(
            sigma_xx, sigma_yy, tau_xy, C[i], nu, thickness, lambda_light, S_i_hat
        )

        synthetic_images[:, :, i, 0] = I0_pol
        synthetic_images[:, :, i, 1] = I45_pol
        synthetic_images[:, :, i, 2] = I90_pol
        synthetic_images[:, :, i, 3] = I135_pol

    return (
        synthetic_images,
        principal_diff,
        theta_p,
        sigma_xx,
        sigma_yy,
        tau_xy,
    )