Source code for ggf.core

import numpy as np

from .matlab_funcs import lscov, legendre
from .sci_funcs import legendrePlm


[docs]def legendre2ggf(coeff, poisson_ratio): """Compute the global geometric factor from Legendre coefficients The definition of the Legendre coefficients is given in :ref:`the theory section <sec_theory_ggf>`. Parameters ---------- coeff: 1d ndarray Legendre coefficients as defined in :cite:`Lure1964` poisson_ratio: float Poisson's ratio of the stretched material. Set this to 0.5 for volume conservation. Returns ------- ggf: float Global geometric factor Notes ----- All odd Legendre coefficients are assumed to be zero, because the stress profile is symmetric with respect to the stretcher axis. """ m = 1 / poisson_ratio def Delta(n): return n * (n - 1) + (2 * n + 1) * (m + 1) / m def L_n(n): return -1 / Delta(n) * (2 * n + 1) * (n + 1) * (n - 2 + 4 / m) def M_n(n): return 1 / Delta(n) * (2 * n + 1) * \ (n**2 + 2 * n - 1 + 2 / m) * n / (n - 1) # Q_n = lambda n: -1/Delta(n) * (2*n+1) * (n + 5 - 4/m) # S_n = lambda n: M_n(n) / n x = 1 # evaluate displacements at the boundary of the sphere theta = 0 # evaluate displacements only on the trapping axis # We use the notation: u_r(theta=0, R=radius) / radius = GGF / G # Thus, in Lur'e eq. (6.6.8), we move radius and G to the left. # ggf = u_r * G / radius ggf = 0 for n, sn in enumerate(coeff): if n == 0: # n=0 contribution: ggf += (m - 2) * sn / (2 * (m + 1)) elif n % 2: if not np.allclose(sn, 0): msg = "Odd coeffecient n={} is non-zero: {}".format(n, sn) raise ValueError(msg) else: ggf += 1 / 8 * 2 * sn / (2 * n + 1) \ * (L_n(n) * x**n + M_n(n) * x**(n - 2)) \ * np.real_if_close(legendre(n, np.cos(theta))[0][0]) # Note that u_theta is not considered here! return ggf
[docs]def stress2legendre(stress, theta, n_poly): """Decompose stress into even Legendre Polynomials The definition of the Legendre decomposition is given in :ref:`the theory section <sec_theory_ggf>`. Parameters ---------- stress: 1d ndarray Radial stress profile (in imaging plane) theta: 1d ndarray Polar angles corresponding to `stress` n_poly: int Number of Legendre polynomials to use Returns ------- coeff: 1d ndarray Legendre coefficients as defined in :cite:`Lure1964` Notes ----- All odd Legendre coefficients are assumed to be zero, because the stress profile is symmetric with respect to the stretcher axis. Therefore, only `n_poly/2` polynomials are considered. """ # Sigma = Sum_n [Coeff(n) P_n(np.cos(theta))] # number of Legendre polynomials used in fit nmax = n_poly # transfer data from stress plot into pair of corresponding variables # [Theta,Sigma] numpoints = theta.shape[0] theta = theta.reshape(-1, 1) sigma = stress.reshape(-1, 1) # Write set of linear equations for stresses in terms of Legendre functions legmat = np.zeros((numpoints, nmax), dtype=float) for ii in range(numpoints): # skip odd Legendre Polynomials since stress is an even function # (symmetrical) for jj in np.arange(nmax)[::2]: legmat[ii, jj] = np.real_if_close( legendrePlm(0, jj, np.cos(theta[ii]))) coeff = lscov(legmat, sigma) return coeff
[docs]def stress2ggf(stress, theta, poisson_ratio, n_poly=120): """Compute the GGf from radial stress using Legendre decomposition Parameters ---------- stress: 1d ndarray Radial stress profile (in imaging plane) theta: 1d ndarray Polar angles corresponding to `stress` poisson_ratio: float Poisson's ratio of the stretched material. Set this to 0.5 for volume conservation. n_poly: int Number of Legendre polynomials to use Returns ------- ggf: float Global geometric factor Notes ----- All odd Legendre coefficients are assumed to be zero, because the stress profile is symmetric with respect to the stretcher axis. Therefore, only `n_poly/2` polynomials are considered. """ coeff = stress2legendre(stress=stress, theta=theta, n_poly=n_poly) ggf = legendre2ggf(coeff=coeff, poisson_ratio=poisson_ratio) return ggf