src.acoustools.BEM.Stiffness

View Source

 1import torch
 2from acoustools.BEM.Force import BEM_compute_force
 3from acoustools.Utilities import create_points, TRANSDUCERS
 4from acoustools.Constants import V
 5import acoustools.Constants as c
 6
 7from torch import Tensor
 8from vedo import Mesh
 9
10def stiffness_finite_differences_BEM(activations:Tensor, points:Tensor, board:Tensor|None=None, scatterer:Mesh = None, path=None, H=None, V=V, delta= 0.001, 
11                                    p_ref=c.P_ref,k=c.k, transducer_radius=c.radius, 
12                                    medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p):
13    '''
14    Computes the stiffness at a point as the gradient of the force. Force computed analytically and then finite differences used to find the gradient \n
15    Computed as `-1* (Fx + Fy + Fz)` where `Fa` is the gradient of force in that direction \n 
16    :param activation: Hologram
17    :param points: Points of interest
18    :param board: Transducers to use
19    :param delta: finite differences step size
20    
21    '''
22
23    if board is None:
24        board = TRANSDUCERS
25
26    dx = create_points(1,1,delta,0,0)
27    dy = create_points(1,1,0,delta,0)
28    dz = create_points(1,1,0,0,delta)
29
30    Fx1 = BEM_compute_force(activations,points + dx,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[0]
31    Fx2 = BEM_compute_force(activations,points - dx,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[0]
32
33    Fx = ((Fx1 - Fx2) / (2*delta))
34
35    Fy1 = BEM_compute_force(activations,points + dy,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[1]
36    Fy2 = BEM_compute_force(activations,points - dy,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[1]
37
38    Fy = ((Fy1 - Fy2) / (2*delta))
39
40    Fz1 = BEM_compute_force(activations,points + dz,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[2]
41    Fz2 = BEM_compute_force(activations,points - dz,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[2]
42    
43    Fz = ((Fz1 - Fz2) / (2*delta))
44
45    return -1* (Fx + Fy + Fz)

def stiffness_finite_differences_BEM( activations: torch.Tensor, points: torch.Tensor, board: torch.Tensor | None = None, scatterer: vedo.mesh.Mesh = None, path=None, H=None, V=4.188790204666667e-09, delta=0.001, p_ref=3.4000000000000004, k=732.7329804081634, transducer_radius=0.0045, medium_density=1.2, medium_speed=343, particle_density=29.36, particle_speed=1052): View Source

11def stiffness_finite_differences_BEM(activations:Tensor, points:Tensor, board:Tensor|None=None, scatterer:Mesh = None, path=None, H=None, V=V, delta= 0.001, 
12                                    p_ref=c.P_ref,k=c.k, transducer_radius=c.radius, 
13                                    medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p):
14    '''
15    Computes the stiffness at a point as the gradient of the force. Force computed analytically and then finite differences used to find the gradient \n
16    Computed as `-1* (Fx + Fy + Fz)` where `Fa` is the gradient of force in that direction \n 
17    :param activation: Hologram
18    :param points: Points of interest
19    :param board: Transducers to use
20    :param delta: finite differences step size
21    
22    '''
23
24    if board is None:
25        board = TRANSDUCERS
26
27    dx = create_points(1,1,delta,0,0)
28    dy = create_points(1,1,0,delta,0)
29    dz = create_points(1,1,0,0,delta)
30
31    Fx1 = BEM_compute_force(activations,points + dx,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[0]
32    Fx2 = BEM_compute_force(activations,points - dx,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[0]
33
34    Fx = ((Fx1 - Fx2) / (2*delta))
35
36    Fy1 = BEM_compute_force(activations,points + dy,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[1]
37    Fy2 = BEM_compute_force(activations,points - dy,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[1]
38
39    Fy = ((Fy1 - Fy2) / (2*delta))
40
41    Fz1 = BEM_compute_force(activations,points + dz,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[2]
42    Fz2 = BEM_compute_force(activations,points - dz,board=board,scatterer=scatterer, path=path, H=H, V=V, p_ref=p_ref, transducer_radius=transducer_radius, medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density, particle_speed=particle_speed)[2]
43    
44    Fz = ((Fz1 - Fz2) / (2*delta))
45
46    return -1* (Fx + Fy + Fz)

Computes the stiffness at a point as the gradient of the force. Force computed analytically and then finite differences used to find the gradient

Computed as -1* (Fx + Fy + Fz) where Fa is the gradient of force in that direction

Parameters

activation: Hologram
points: Points of interest
board: Transducers to use
delta: finite differences step size