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