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, propagate
  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_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:float=c.V, scatterer:Mesh=None, 
 17                  H:Tensor=None, path:str="Media", k=c.k, transducer_radius=c.radius, p_ref=c.P_ref,
 18                  medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p) -> Tensor | tuple[Tensor, Tensor, Tensor]:
 19    '''
 20    Returns the force on a particle using the analytical derivative of the Gor'kov potential and BEM\n
 21    :param activations: Transducer hologram
 22    :param points: Points to propagate to
 23    :param board: Transducers to use, if `None` uses `acoustools.Utilities.TRANSDUCERS`
 24    :param return_components: If true returns force as one tensor otherwise returns Fx, Fy, Fz
 25    :param V: Particle volume
 26    :param scatterer: Scatterer to use
 27    :param H: H to use, will load/compute if None
 28    :param path: Path to folder containing BEMCache
 29    :return: force  
 30    '''
 31
 32    #Bk.2 Pg.319
 33
 34    if board is None:
 35        board = TRANSDUCERS
 36    
 37    F = compute_E(scatterer=scatterer,points=points,board=board,H=H, path=path, k=k, transducer_radius=transducer_radius, p_ref=p_ref)
 38    Fx, Fy, Fz = BEM_forward_model_grad(points,transducers=board,scatterer=scatterer,H=H, path=path, k=k, transducer_radius=transducer_radius, p_ref=p_ref)
 39    Fxx, Fyy, Fzz = BEM_forward_model_second_derivative_unmixed(points,transducers=board,scatterer=scatterer,H=H, path=path, k=k, transducer_radius=transducer_radius, p_ref=p_ref)
 40    Fxy, Fxz, Fyz = BEM_forward_model_second_derivative_mixed(points,transducers=board,scatterer=scatterer,H=H, path=path, k=k, transducer_radius=transducer_radius, p_ref=p_ref)
 41
 42    p   = (F@activations)
 43    Px  = (Fx@activations)
 44    Py  = (Fy@activations)
 45    Pz  = (Fz@activations)
 46    Pxx = (Fxx@activations)
 47    Pyy = (Fyy@activations)
 48    Pzz = (Fzz@activations)
 49    Pxy = (Fxy@activations)
 50    Pxz = (Fxz@activations)
 51    Pyz = (Fyz@activations)
 52
 53
 54    grad_p = torch.stack([Px,Py,Pz])
 55    grad_px = torch.stack([Pxx,Pxy,Pxz])
 56    grad_py = torch.stack([Pxy,Pyy,Pyz])
 57    grad_pz = torch.stack([Pxz,Pyz,Pzz])
 58
 59
 60    p_term = p*grad_p.conj() + p.conj()*grad_p
 61
 62    px_term = Px*grad_px.conj() + Px.conj()*grad_px
 63    py_term = Py*grad_py.conj() + Py.conj()*grad_py
 64    pz_term = Pz*grad_pz.conj() + Pz.conj()*grad_pz
 65
 66    K1, K2 = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
 67    
 68    grad_U = K1 * p_term - K2 * (px_term + py_term + pz_term)
 69    force = -1*(grad_U).squeeze().real
 70
 71    if return_components:
 72        return force[0], force[1], force[2] 
 73    else:
 74        return force 
 75
 76def torque_mesh_surface(activations:Tensor, scatterer:Mesh=None, board:Tensor|None=None, sum_elements = True, use_pressure:bool=False,
 77                       H:Tensor=None, diameter=c.wavelength*2, p_ref = c.P_ref,
 78                       path:str="Media", surface_path:str = "/Sphere-solidworks-lam2.stl",
 79                       surface:Mesh|None=None, use_cache_H:bool=True, 
 80                       E:Tensor|None=None, Ex:Tensor|None=None, Ey:Tensor|None=None, Ez:Tensor|None=None,
 81                        internal_points=None ) -> Tensor | tuple[Tensor, Tensor, Tensor]:
 82    '''
 83    Computes the force on a scattering obejct by computing thr force on a far field surface\\
 84    :param activations: Hologram
 85    :param scatterer: Object to compute force on
 86    :param board: Transducers
 87    :param sum_elements: if True will call sum across mesh elements
 88    :param H: H matrix to use for BEM, if None will be computed
 89    :param diameter: diameter of surfac to use
 90    :param path: path to BEMCache
 91    :param surface_path: Name of stl to use for surface such that path + surface path is the full address
 92    :param surface: Surface to use, if None will be laoded from surface_path
 93    :param use_cache_H: If true use BEM cache system
 94    :param E: E matrix to use for BEM, if None will be computed
 95    :param Ex: Grad E wrt x matrix to use for BEM, if None will be computed
 96    :param Ey: Grad E wrt y matrix to use for BEM, if None will be computed
 97    :param Ey: Grad E wrt z matrix to use for BEM, if None will be computed
 98    :returns torque:
 99    '''
100
101    object_com = get_centre_of_mass_as_points(scatterer)
102
103    if surface is None:
104        surface = load_scatterer(path+surface_path)
105        scale_to_diameter(surface,diameter, reset=False, origin=True)
106        centre_scatterer(surface)
107        translate(surface, dx = object_com[:,0].item(), dy=object_com[:,1].item(), dz = object_com[:,2].item())
108    
109
110
111    points = get_centres_as_points(surface)
112    norms = get_normals_as_points(surface)
113    areas = get_areas(surface)
114    surface_com = get_centre_of_mass_as_points(surface)
115    r = points - surface_com
116    # print('dot',torch.sum(r * norms, dim=1) / torch.norm(r,2, dim=1)) 
117    
118    if use_pressure:
119        if E is None:
120            E = compute_E(scatterer, points, board, use_cache_H=use_cache_H, path=path, H=H,internal_points=internal_points, p_ref=p_ref)
121        p = propagate(activations,points,board,A=E, p_ref=p_ref)
122        pressure_square = torch.abs(p)**2
123        pressure_time_average = 1/2 * pressure_square
124
125        r_cross_n = torch.cross(r,norms, dim=1)
126        pressure_term = pressure_time_average * r_cross_n * areas
127        if sum_elements: pressure_term = torch.sum(pressure_term, dim=2)
128    else:
129        pressure_term = 0
130
131    if Ex is None or Ey is None or Ez is None:
132        Ex, Ey, Ez = BEM_forward_model_grad(points, scatterer, board, use_cache_H=use_cache_H, H=H, path=path,internal_points=internal_points, p_ref=p_ref)
133    
134    px = (Ex@activations).squeeze(2).unsqueeze(0)
135    py = (Ey@activations).squeeze(2).unsqueeze(0)
136    pz = (Ez@activations).squeeze(2).unsqueeze(0)
137
138    grad  = torch.cat((px,py,pz),dim=1)
139    velocity = grad /( 1j * c.p_0 * c.angular_frequency)
140    # return torch.sum(velocity,dim=2)
141    # r_cross_outer_conj = torch.cross(r, outer_conj)
142
143    # time_average_outer_v = 1/2 * torch.einsum('bin,bjn -> bijn', r_cross_v, velocity.conj().resolve_conj()) #Takes two (B,3,N) vectors and computes the outer product on them - i think...
144    vconj_dot_n = torch.sum(velocity.conj() * norms, dim=1, keepdim=True)
145    time_average_velocity = 1/2 *(torch.cross((vconj_dot_n * r), velocity, dim=1)).real * areas
146    if sum_elements: time_average_velocity = torch.sum(time_average_velocity, dim=2)
147    torque = - pressure_term - c.p_0 * time_average_velocity
148
149    return torque
150
151
152
153def force_mesh_surface(activations:Tensor, scatterer:Mesh=None, board:Tensor|None=None,
154                       return_components:bool=False, sum_elements:bool = True, return_momentum:bool = False,
155                       H:Tensor=None, diameter:float=c.wavelength*2,
156                       path:str="Media", surface_path:str = "/Sphere-solidworks-lam2.stl",
157                       surface:Mesh|None=None, use_cache_H:bool=True, 
158                       E:Tensor|None=None, Ex:Tensor|None=None, Ey:Tensor|None=None, Ez:Tensor|None=None, 
159                       use_momentum:float=True, p_ref:float=c.P_ref, internal_points=None, k=c.k, transducer_radius=c.radius) -> Tensor | tuple[Tensor, Tensor, Tensor]:
160    '''
161    Computes the torque on a scattering obejct by computing thr force on a far field surface\\
162    :param activations: Hologram
163    :param scatterer: Object to compute force on
164    :param board: Transducers
165    :param return_components: if True will return Fx,Fy,Fz else returns force
166    :param sum_elements: if True will call sum across mesh elements
167    :param return_momentum: if True will return total force and the momentum flux term
168    :param H: H matrix to use for BEM, if None will be computed
169    :param diameter: diameter of surface to use
170    :param path: path to BEMCache
171    :param surface_path: Name of stl to use for surface such that path + surface path is the full address
172    :param surface: Surface to use, if None will be laoded from surface_path
173    :param use_cache_H: If true use BEM cache system
174    :param E: E matrix to use for BEM, if None will be computed
175    :param Ex: Grad E wrt x matrix to use for BEM, if None will be computed
176    :param Ey: Grad E wrt y matrix to use for BEM, if None will be computed
177    :param Ey: Grad E wrt z matrix to use for BEM, if None will be computed
178    :use_momentum: If true will use momentum flux terms
179    :returns force:
180    '''
181    
182    if surface is None:
183        surface = load_scatterer(path+surface_path)
184        scale_to_diameter(surface,diameter, reset=False, origin=True)
185        centre_scatterer(surface)
186        object_com = get_centre_of_mass_as_points(scatterer)
187        translate(surface, dx = object_com[:,0].item(), dy=object_com[:,1].item(), dz = object_com[:,2].item())
188        
189
190    points = get_centres_as_points(surface)
191    norms = get_normals_as_points(surface)
192    areas = get_areas(surface)
193    
194    if E is None:
195        E,F,G,H = compute_E(scatterer, points, board,path=path, H=H, return_components=True, use_cache_H=use_cache_H, p_ref=p_ref,internal_points=internal_points, k=k, transducer_radius=transducer_radius)
196    
197    force, momentum = force_mesh(activations, points,norms,areas,board=board,F=E, use_momentum=use_momentum, p_ref=p_ref, k=k, transducer_radius=transducer_radius,
198                    grad_function=BEM_forward_model_grad, grad_function_args={'scatterer':scatterer,
199                                                                                'H':H,
200                                                                                'path':path,
201                                                                                "internal_points":internal_points,
202                                                                                }, 
203                                                                                return_components=True,
204                                                                                Ax = Ex, Ay=Ey, Az=Ez)
205
206
207
208    if sum_elements: force=torch.sum(force, dim=2)
209
210    if return_components:
211        if not return_momentum:
212            return (force[:,0]), (force[:,1]), (force[:,2])
213        else:
214            return (force[:,0]), (force[:,1]), (force[:,2]), momentum
215   
216    if return_momentum: return force, momentum
217    return force
218
219def force_mesh_surface_divergance(activations:Tensor, scatterer:Mesh=None, board:Tensor|None=None,
220                       sum_elements:bool = True, H:Tensor=None, diameter:float=c.wavelength*2,
221                       path:str="Media", surface_path:str = "/Sphere-solidworks-lam2.stl",
222                       surface:Mesh|None=None, use_cache_H:bool=True, force:Tensor|None=None, norms:Tensor|None=None, p_ref=c.P_ref, k=c.k, transducer_radius=c.radius) -> Tensor:
223    '''
224    Computes the divergance force (the dot product of the force and normals on the surface) on a scattering obejct by computing thr force on a far field surface\\
225    :param activations: Hologram
226    :param scatterer: Object to compute force on
227    :param board: Transducers
228    :param sum_elements: if True will call sum across mesh elements
229    :param H: H matrix to use for BEM, if None will be computed
230    :param diameter: diameter of surfac to use
231    :param path: path to BEMCache
232    :param surface_path: Name of stl to use for surface such that path + surface path is the full address
233    :param surface: Surface to use, if None will be laoded from surface_path
234    :param use_cache_H: If true use BEM cache system
235    :param force: Precomputed force, if non is computed using `force_mesh_surface`
236    :param norms: Precomputed norms, if non is found from surface
237    :returns divergance of force:
238    '''
239
240
241    if surface is None:
242        surface = load_scatterer(path+surface_path)
243        scale_to_diameter(surface,diameter, reset=False, origin=True)
244        centre_scatterer(surface)
245        object_com = get_centre_of_mass_as_points(scatterer)
246        translate(surface, dx = object_com[:,0].item(), dy=object_com[:,1].item(), dz = object_com[:,2].item())
247    
248    if force is None:
249        force = force_mesh_surface(activations, scatterer, board, H=H, diameter=diameter, path=path, 
250                               surface_path=surface_path, surface=surface, use_cache_H=use_cache_H, sum_elements=False, use_momentum=True, k=k, p_ref=p_ref, transducer_radius=transducer_radius) 
251
252    if norms is None: norms = get_normals_as_points(surface)
253    areas = get_areas(surface)
254
255    div = (torch.sum(norms * force, dim=1) * areas )
256
257    if sum_elements: div = torch.sum(div, dim=1)
258
259    v = get_volume(surface)
260
261    return div / v
262
263
264def force_mesh_surface_curl(activations:Tensor, scatterer:Mesh=None, board:Tensor|None=None,
265                       sum_elements:bool = True, H:Tensor=None, diameter:float=c.wavelength*2,
266                       path:str="Media", surface_path:str = "/Sphere-solidworks-lam2.stl",
267                       surface:Mesh|None=None, use_cache_H:bool=True, magnitude:Tensor|None = False, force:Tensor|None=None, p_ref=c.P_ref, k=c.k, transducer_radius=c.radius) -> Tensor:
268    '''
269    Computes the curl force (the cross product of the force and normals on the surface) on a scattering obejct by computing thr force on a far field surface\\
270    :param activations: Hologram
271    :param scatterer: Object to compute force on
272    :param board: Transducers
273    :param sum_elements: if True will call sum across mesh elements
274    :param H: H matrix to use for BEM, if None will be computed
275    :param diameter: diameter of surfac to use
276    :param path: path to BEMCache
277    :param surface_path: Name of stl to use for surface such that path + surface path is the full address
278    :param surface: Surface to use, if None will be laoded from surface_path
279    :param use_cache_H: If true use BEM cache system
280    :param force: Precomputed force, if non is computed using `force_mesh_surface`
281    :param magnitude: If true will call `torch.norm` on the curl vectors
282    :returns curl of force:
283    '''
284    
285    if force is None: force = force_mesh_surface(activations, scatterer, board, H=H, diameter=diameter, path=path, 
286                               surface_path=surface_path, surface=surface, use_cache_H=use_cache_H, sum_elements=False, k=k, p_ref=p_ref, transducer_radius=transducer_radius) 
287
288    if surface is None:
289        surface = load_scatterer(path+surface_path)
290        scale_to_diameter(surface,diameter, reset=False, origin=True)
291        centre_scatterer(surface)
292        object_com = get_centre_of_mass_as_points(scatterer)
293        translate(surface, dx = object_com[:,0].item(), dy=object_com[:,1].item(), dz = object_com[:,2].item())
294        
295    norms = get_normals_as_points(surface).real
296    areas = get_areas(surface)
297
298    curl = torch.cross(force, norms, dim=1) * areas 
299
300    if sum_elements: curl = torch.sum(curl, dim=2)
301
302    v = get_volume(surface)
303
304    curl = curl/v
305
306    if magnitude: return torch.norm(curl, dim=1)
307
308    return curl 
309
310
311def get_force_mesh_along_axis(start:Tensor,end:Tensor, activations:Tensor, scatterers:list[Mesh], board:Tensor, mask:Tensor|None=None, steps:int=200, 
312                              path:str="Media",print_lines:bool=False, use_cache:bool=True, 
313                              Hs:Tensor|None = None, Hxs:Tensor|None=None, Hys:Tensor|None=None, Hzs:Tensor|None=None, p_ref=c.P_ref, k=c.k, transducer_radius=c.radius) -> tuple[list[Tensor],list[Tensor],list[Tensor]]:
314    '''
315    @private
316    Computes the force on a mesh at each point from `start` to `end` with number of samples = `steps`  \n
317    :param start: The starting position
318    :param end: The ending position
319    :param activations: Transducer hologram
320    :param scatterers: First element is the mesh to move, rest is considered static reflectors 
321    :param board: Transducers to use 
322    :param mask: The mask to apply to filter force for only the mesh to move
323    :param steps: Number of steps to take from start to end
324    :param path: path to folder containing BEMCache/ 
325    :param print_lines: if true prints messages detaling progress
326    :param use_cache: If true uses the cache system, otherwise computes H and does not save it
327    :param Hs: List of precomputed forward propagation matricies
328    :param Hxs: List of precomputed derivative of forward propagation matricies wrt x
329    :param Hys: List of precomputed derivative of forward propagation matricies wrt y
330    :param Hzs: List of precomputed derivative of forward propagation matricies wrt z
331    :return: list for each axis of the force at each position
332    '''
333
334    print("get_force_mesh_along_axis - implementation H grad incorrect - do not use")
335    # if Ax is None or Ay is None or Az is None:
336    #     Ax, Ay, Az = grad_function(points=points, transducers=board, **grad_function_args)
337    direction = (end - start) / steps
338
339    translate(scatterers[0], start[0].item() - direction[0].item(), start[1].item() - direction[1].item(), start[2].item() - direction[2].item())
340    scatterer = merge_scatterers(*scatterers)
341
342    points = get_centres_as_points(scatterer)
343    if mask is None:
344        mask = torch.ones(points.shape[2]).to(bool)
345
346    Fxs = []
347    Fys = []
348    Fzs = []
349
350    for i in range(steps+1):
351        if print_lines:
352            print(i)
353        
354        
355        translate(scatterers[0], direction[0].item(), direction[1].item(), direction[2].item())
356        scatterer = merge_scatterers(*scatterers)
357
358        points = get_centres_as_points(scatterer)
359        areas = get_areas(scatterer)
360        norms = get_normals_as_points(scatterer)
361
362        if Hs is None:
363            H = get_cache_or_compute_H(scatterer, board, path=path, print_lines=print_lines, use_cache_H=use_cache, k=k, p_ref=p_ref, transducer_radius=transducer_radius)
364        else:
365            H = Hs[i]
366        
367        if Hxs is None or Hys is None or Hzs is None:
368            Hx, Hy, Hz = get_cache_or_compute_H_gradients(scatterer, board, path=path, print_lines=print_lines, use_cache_H_grad=use_cache)
369        else:
370            Hx = Hxs[i]
371            Hy = Hys[i]
372            Hz = Hzs[i]
373        
374
375        force = force_mesh(activations, points, norms, areas, board, F=H, Ax=Hx, Ay=Hy, Az=Hz)
376
377        force = torch.sum(force[:,:,mask],dim=2).squeeze()
378        Fxs.append(force[0])
379        Fys.append(force[1])
380        Fzs.append(force[2])
381        
382        # print(i, force[0].item(), force[1].item(),force[2].item())
383    return Fxs, Fys, Fzs