src.acoustools.BEM.Force

  1from acoustools.BEM import BEM_forward_model_second_derivative_mixed, BEM_forward_model_second_derivative_unmixed, BEM_forward_model_grad, compute_E, get_cache_or_compute_H, get_cache_or_compute_H_gradients
  2from acoustools.Utilities import TRANSDUCERS
  3from acoustools.Force import force_mesh
  4from acoustools.Mesh import load_scatterer, get_centres_as_points, get_normals_as_points, get_areas, scale_to_diameter,\
  5    centre_scatterer, translate, merge_scatterers, get_edge_data, get_centre_of_mass_as_points, get_volume
  6from acoustools.Gorkov import get_gorkov_constants
  7
  8
  9import acoustools.Constants as c
 10
 11from torch import Tensor
 12import torch
 13
 14from vedo import Mesh
 15
 16def BEM_compute_force(activations:Tensor, points:Tensor,board:Tensor|None=None,return_components:bool=False, V=c.V, scatterer:Mesh=None, 
 17                  H:Tensor=None, path:str="Media") -> Tensor | tuple[Tensor, Tensor, Tensor]:
 18    '''
 19    Returns the force on a particle using the analytical derivative of the Gor'kov potential and BEM\n
 20    :param activations: Transducer hologram
 21    :param points: Points to propagate to
 22    :param board: Transducers to use, if `None` uses `acoustools.Utilities.TRANSDUCERS`
 23    :param return_components: If true returns force as one tensor otherwise returns Fx, Fy, Fz
 24    :param V: Particle volume
 25    :param scatterer: Scatterer to use
 26    :param H: H to use, will load/compute if None
 27    :param path: Path to folder containing BEMCache
 28    :return: force  
 29    '''
 30
 31    #Bk.2 Pg.319
 32
 33    if board is None:
 34        board = TRANSDUCERS
 35    
 36    F = compute_E(scatterer=scatterer,points=points,board=board,H=H, path=path)
 37    Fx, Fy, Fz = BEM_forward_model_grad(points,transducers=board,scatterer=scatterer,H=H, path=path)
 38    Fxx, Fyy, Fzz = BEM_forward_model_second_derivative_unmixed(points,transducers=board,scatterer=scatterer,H=H, path=path)
 39    Fxy, Fxz, Fyz = BEM_forward_model_second_derivative_mixed(points,transducers=board,scatterer=scatterer,H=H, path=path)
 40
 41    p   = (F@activations)
 42    Px  = (Fx@activations)
 43    Py  = (Fy@activations)
 44    Pz  = (Fz@activations)
 45    Pxx = (Fxx@activations)
 46    Pyy = (Fyy@activations)
 47    Pzz = (Fzz@activations)
 48    Pxy = (Fxy@activations)
 49    Pxz = (Fxz@activations)
 50    Pyz = (Fyz@activations)
 51
 52
 53    grad_p = torch.stack([Px,Py,Pz])
 54    grad_px = torch.stack([Pxx,Pxy,Pxz])
 55    grad_py = torch.stack([Pxy,Pyy,Pyz])
 56    grad_pz = torch.stack([Pxz,Pyz,Pzz])
 57
 58
 59    p_term = p*grad_p.conj() + p.conj()*grad_p
 60
 61    px_term = Px*grad_px.conj() + Px.conj()*grad_px
 62    py_term = Py*grad_py.conj() + Py.conj()*grad_py
 63    pz_term = Pz*grad_pz.conj() + Pz.conj()*grad_pz
 64
 65    K1, K2 = get_gorkov_constants(V=V)
 66    
 67    grad_U = K1 * p_term - K2 * (px_term + py_term + pz_term)
 68    force = -1*(grad_U).squeeze().real
 69
 70    if return_components:
 71        return force[0], force[1], force[2] 
 72    else:
 73        return force 
 74
 75
 76def force_mesh_surface(activations:Tensor, scatterer:Mesh=None, board:Tensor|None=None,
 77                       return_components:bool=False, sum_elements = True, return_momentum = False,
 78                       H:Tensor=None, diameter=c.wavelength*2,
 79                       path:str="Media", surface_path:str = "/Sphere-solidworks-lam2.stl",
 80                       surface:Mesh|None=None, use_cache_H:bool=True, E=None, Ex=None, Ey=None, Ez=None, use_momentum=True) -> Tensor | tuple[Tensor, Tensor, Tensor]:
 81    
 82    if surface is None:
 83        surface = load_scatterer(path+surface_path)
 84        scale_to_diameter(surface,diameter, reset=False, origin=True)
 85        centre_scatterer(surface)
 86        object_com = get_centre_of_mass_as_points(scatterer)
 87        translate(surface, dx = object_com[:,0].item(), dy=object_com[:,1].item(), dz = object_com[:,2].item())
 88        
 89
 90    points = get_centres_as_points(surface)
 91    norms = get_normals_as_points(surface)
 92    areas = get_areas(surface)
 93    
 94    if E is None:
 95        E,F,G,H = compute_E(scatterer, points, board,path=path, H=H, return_components=True, use_cache_H=use_cache_H)
 96    
 97    force, momentum = force_mesh(activations, points,norms,areas,board=board,F=E, use_momentum=use_momentum,
 98                    grad_function=BEM_forward_model_grad, grad_function_args={'scatterer':scatterer,
 99                                                                                'H':H,
100                                                                                'path':path}, return_components=True,
101                                                                                Ax = Ex, Ay=Ey, Az=Ez)
102
103
104
105    if sum_elements: force=torch.sum(force, dim=2)
106
107    if return_components:
108        if not return_momentum:
109            return (force[:,0]), (force[:,1]), (force[:,2])
110        else:
111            return (force[:,0]), (force[:,1]), (force[:,2]), momentum
112   
113    if return_momentum: return force, momentum
114    return force
115
116def force_mesh_surface_divergance(activations:Tensor, scatterer:Mesh=None, board:Tensor|None=None,
117                       sum_elements = True, H:Tensor=None, diameter=c.wavelength*2,
118                       path:str="Media", surface_path:str = "/Sphere-solidworks-lam2.stl",
119                       surface:Mesh|None=None, use_cache_H:bool=True, force=None, norms=None) -> Tensor:
120
121
122    if surface is None:
123        surface = load_scatterer(path+surface_path)
124        scale_to_diameter(surface,diameter, reset=False, origin=True)
125        centre_scatterer(surface)
126        object_com = get_centre_of_mass_as_points(scatterer)
127        translate(surface, dx = object_com[:,0].item(), dy=object_com[:,1].item(), dz = object_com[:,2].item())
128    
129    if force is None:
130        force = force_mesh_surface(activations, scatterer, board, H=H, diameter=diameter, path=path, 
131                               surface_path=surface_path, surface=surface, use_cache_H=use_cache_H, sum_elements=False, use_momentum=True) 
132
133    if norms is None: norms = get_normals_as_points(surface)
134    areas = get_areas(surface)
135
136    div = (torch.sum(norms * force, dim=1) * areas )
137
138    if sum_elements: div = torch.sum(div, dim=1)
139
140    v = get_volume(surface)
141
142    return div / v
143
144
145def force_mesh_surface_curl(activations:Tensor, scatterer:Mesh=None, board:Tensor|None=None,
146                       sum_elements = True, H:Tensor=None, diameter=c.wavelength*2,
147                       path:str="Media", surface_path:str = "/Sphere-solidworks-lam2.stl",
148                       surface:Mesh|None=None, use_cache_H:bool=True, magnitude = False, force=None) -> Tensor:
149    
150    if force is None: force = force_mesh_surface(activations, scatterer, board, H=H, diameter=diameter, path=path, 
151                               surface_path=surface_path, surface=surface, use_cache_H=use_cache_H, sum_elements=False) 
152
153    if surface is None:
154        surface = load_scatterer(path+surface_path)
155        scale_to_diameter(surface,diameter, reset=False, origin=True)
156        centre_scatterer(surface)
157        object_com = get_centre_of_mass_as_points(scatterer)
158        translate(surface, dx = object_com[:,0].item(), dy=object_com[:,1].item(), dz = object_com[:,2].item())
159        
160    norms = get_normals_as_points(surface).real
161    areas = get_areas(surface)
162
163    curl = torch.cross(force, norms, dim=1) * areas 
164
165    if sum_elements: curl = torch.sum(curl, dim=2)
166
167    v = get_volume(surface)
168
169    curl = curl/v
170
171    if magnitude: return torch.norm(curl, dim=1)
172
173    return curl 
174
175
176def get_force_mesh_along_axis(start:Tensor,end:Tensor, activations:Tensor, scatterers:list[Mesh], board:Tensor, mask:Tensor|None=None, steps:int=200, 
177                              path:str="Media",print_lines:bool=False, use_cache:bool=True, 
178                              Hs:Tensor|None = None, Hxs:Tensor|None=None, Hys:Tensor|None=None, Hzs:Tensor|None=None) -> tuple[list[Tensor],list[Tensor],list[Tensor]]:
179    '''
180    Computes the force on a mesh at each point from `start` to `end` with number of samples = `steps`  \n
181    :param start: The starting position
182    :param end: The ending position
183    :param activations: Transducer hologram
184    :param scatterers: First element is the mesh to move, rest is considered static reflectors 
185    :param board: Transducers to use 
186    :param mask: The mask to apply to filter force for only the mesh to move
187    :param steps: Number of steps to take from start to end
188    :param path: path to folder containing BEMCache/ 
189    :param print_lines: if true prints messages detaling progress
190    :param use_cache: If true uses the cache system, otherwise computes H and does not save it
191    :param Hs: List of precomputed forward propagation matricies
192    :param Hxs: List of precomputed derivative of forward propagation matricies wrt x
193    :param Hys: List of precomputed derivative of forward propagation matricies wrt y
194    :param Hzs: List of precomputed derivative of forward propagation matricies wrt z
195    :return: list for each axis of the force at each position
196    '''
197    # if Ax is None or Ay is None or Az is None:
198    #     Ax, Ay, Az = grad_function(points=points, transducers=board, **grad_function_args)
199    direction = (end - start) / steps
200
201    translate(scatterers[0], start[0].item() - direction[0].item(), start[1].item() - direction[1].item(), start[2].item() - direction[2].item())
202    scatterer = merge_scatterers(*scatterers)
203
204    points = get_centres_as_points(scatterer)
205    if mask is None:
206        mask = torch.ones(points.shape[2]).to(bool)
207
208    Fxs = []
209    Fys = []
210    Fzs = []
211
212    for i in range(steps+1):
213        if print_lines:
214            print(i)
215        
216        
217        translate(scatterers[0], direction[0].item(), direction[1].item(), direction[2].item())
218        scatterer = merge_scatterers(*scatterers)
219
220        points = get_centres_as_points(scatterer)
221        areas = get_areas(scatterer)
222        norms = get_normals_as_points(scatterer)
223
224        if Hs is None:
225            H = get_cache_or_compute_H(scatterer, board, path=path, print_lines=print_lines, use_cache_H=use_cache)
226        else:
227            H = Hs[i]
228        
229        if Hxs is None or Hys is None or Hzs is None:
230            Hx, Hy, Hz = get_cache_or_compute_H_gradients(scatterer, board, path=path, print_lines=print_lines, use_cache_H_grad=use_cache)
231        else:
232            Hx = Hxs[i]
233            Hy = Hys[i]
234            Hz = Hzs[i]
235        
236
237        force = force_mesh(activations, points, norms, areas, board, F=H, Ax=Hx, Ay=Hy, Az=Hz)
238
239        force = torch.sum(force[:,:,mask],dim=2).squeeze()
240        Fxs.append(force[0])
241        Fys.append(force[1])
242        Fzs.append(force[2])
243        
244        # print(i, force[0].item(), force[1].item(),force[2].item())
245    return Fxs, Fys, Fzs