src.acoustools.BEM.Gradients.E_Gradient

View Source

  1import torch
  2from torch import Tensor
  3
  4from vedo import Mesh
  5
  6import hashlib, pickle
  7
  8from acoustools.Utilities import DTYPE, forward_model_grad, forward_model_second_derivative_unmixed, forward_model_second_derivative_mixed
  9from acoustools.BEM.Forward_models import get_cache_or_compute_H
 10from acoustools.Mesh import get_areas, get_centres_as_points, get_normals_as_points
 11import acoustools.Constants as Constants
 12
 13def BEM_forward_model_grad(points:Tensor, scatterer:Mesh, transducers:Tensor=None, use_cache_H:bool=True, 
 14                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False,k=Constants.k,
 15                           path:str="Media", p_ref=Constants.P_ref, internal_points = None, transducer_radius=Constants.radius) -> tuple[Tensor, Tensor, Tensor] | tuple[Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor]:
 16    '''
 17    Computes the gradient of the forward propagation for BEM\n
 18    :param scatterer: The mesh used (as a `vedo` `mesh` object)
 19    :param transducers: Transducers to use, if `None` uses `acoustools.Utilities.TRANSDUCERS`
 20    :param use_cache_H_grad: If true uses the cache system, otherwise computes `H` and does not save it
 21    :param print_lines: if true prints messages detaling progress
 22    :param H: Precomputed `H` - if `None` `H` will be computed
 23    :param return_components: if true will return the subparts used to compute
 24    :param path: path to folder containing `BEMCache/` 
 25    :return: Ex, Ey, Ez
 26    '''
 27    if transducers is None:
 28        transducers = TRANSDUCERS
 29
 30    B = points.shape[0]
 31    if H is None:
 32        H = get_cache_or_compute_H(scatterer,transducers,use_cache_H, path, print_lines,p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius)
 33    
 34    Fx, Fy, Fz  = forward_model_grad(points, transducers,p_ref=p_ref, k=k, transducer_radius=transducer_radius)
 35    Gx, Gy, Gz = get_G_partial(points, scatterer, transducers, k=k)
 36
 37
 38    Ex = Fx + Gx@H 
 39    Ey = Fy + Gy@H
 40    Ez = Fz + Gz@H
 41
 42
 43    if return_components:
 44        return Ex.to(DTYPE), Ey.to(DTYPE), Ez.to(DTYPE), Fx, Fy, Fz, Gx, Gy, Gz, H
 45    else:
 46        return Ex.to(DTYPE), Ey.to(DTYPE), Ez.to(DTYPE)
 47    
 48
 49def get_G_second_unmixed(points:Tensor, scatterer:Mesh, board:Tensor|None=None, return_components:bool=False, k=Constants.k) -> tuple[Tensor, Tensor, Tensor]:
 50
 51    #Bk.3 Pg.80
 52    #P is much higher than the others?
 53
 54    areas = get_areas(scatterer)
 55    centres = get_centres_as_points(scatterer)
 56    normals = get_normals_as_points(scatterer).unsqueeze(2)
 57
 58    N = points.shape[2]
 59    M = centres.shape[2]
 60
 61    points.requires_grad_()
 62    grad_Gx, grad_Gy, grad_Gz, G,P,C,Ga,Pa,Ca = get_G_partial(points, scatterer, board, True, k=k)
 63    # exit()
 64
 65    points  = points.unsqueeze(3)  # [B, 3, N, 1]
 66    centres = centres.unsqueeze(2)  # [B, 3, 1, M]
 67    
 68
 69    diff = points - centres    
 70    diff_square = diff**2
 71    diff_square_sum = torch.sum(diff_square, 1)
 72    distances = torch.sqrt(diff_square_sum)
 73    distances_cube = distances ** 3
 74    distances_five = distances ** 5
 75
 76    distances_expanded = distances.unsqueeze(1).expand((1,3,N,M))
 77    distances_expanded_square = distances_expanded**2
 78    distances_expanded_cube = distances_expanded**3
 79    distances_expanded_four = distances_expanded**4
 80    distances_expanded_five = distances_expanded**5
 81
 82    kd = k * distances_expanded
 83
 84    exp_ikd = torch.exp(1j*kd)
 85    
 86    Gaa = areas/(4*Constants.pi * 1j) *exp_ikd * ((1j*k**2 * diff_square)/ distances_expanded_cube + (1j*k)/distances_expanded_square - 1/distances_expanded_cube - (3*1j*k*diff_square)/distances_expanded_four + (3*diff_square)/distances_expanded_five)
 87    Gxx = Gaa[:,0]
 88    Gyy = Gaa[:,1]
 89    Gzz = Gaa[:,2]
 90    
 91    #P = ik - 1/d
 92    # Paa = (distances_expanded * daa - 2*da**2) / distances_expanded_cube
 93    # Paa = (3 * diff_square - distances_expanded_square) / distances_expanded_five
 94    Paa = 1/distances_expanded_cube - (3*diff_square) / distances_expanded_five  #Equal to lines above
 95    Pxx = Paa[:,0]
 96    Pyy = Paa[:,1]
 97    Pzz = Paa[:,2]
 98
 99    #C = (diff dot c )/ d
100
101    dx = diff[:,0,:,:]
102    dy = diff[:,1,:,:]
103    dz = diff[:,2,:,:]
104
105
106    nx = normals[:,0,:]
107    ny = normals[:,1,:]
108    nz = normals[:,2,:]
109
110    nd = nx * dx + ny * dy + nz * dz
111
112    Cxx = (2 * nx * dx)/distances_cube - nd/distances_cube + (3*nd*dx**2)/distances_five
113    Cyy = (2 * ny * dy)/distances_cube - nd/distances_cube + (3*nd*dy**2)/distances_five
114    Czz = (2 * nz * dz)/distances_cube - nd/distances_cube + (3*nd*dz**2)/distances_five
115
116    # Caa = torch.stack([Cxx, Cyy, Czz], dim=1)
117
118    # Cxx = (((3 * dax**2)/distances_five) - (1/distances_cube)) * nd - (2*nx*dx) / distances_cube
119    # Cyy = (((3 * day**2)/distances_five) - (1/distances_cube)) * nd - (2*ny*dy) / distances_cube 
120    # Czz = (((3 * daz**2)/distances_five) - (1/distances_cube)) * nd - (2*nz*dz) / distances_cube 
121    # Caa = torch.stack([Cxx, Cyy, Czz], dim=1) 
122
123    # grad_2_G_unmixed_old = 2 * Ca * (Ga * P + G * Pa) + C * (2 * Ga*Pa + Gaa * P + G * Paa) + G*P*Caa
124    # grad_2_G_unmixed = C * (2*Ga * Pa + Gaa * P + G*Paa) + P * (2*Ga * Ca + G*Caa) + 2*G * Pa * Ca
125
126    Cx = Ca[:,0]
127    Cy = Ca[:,1]
128    Cz = Ca[:,2]
129
130    Px = Pa[:,0]
131    Py = Pa[:,1]
132    Pz = Pa[:,2]
133
134    Gx = Ga[:,0]
135    Gy = Ga[:,1]
136    Gz = Ga[:,2]
137
138    G = G[:,0,:]
139    P = P[:,0,:] 
140    
141
142
143    grad_2_G_unmixed_xx = Gxx * P * C + G*Pxx*C + G*P*Cxx+  2*Gx*Px*C + 2*Gx*P*Cx + 2*G*Px*Cx
144    grad_2_G_unmixed_yy = Gyy * P * C + G*Pyy*C + G*P*Cyy+  2*Gy*Py*C + 2*Gy*P*Cy + 2*G*Py*Cy
145    grad_2_G_unmixed_zz = Gzz * P * C + G*Pzz*C + G*P*Czz+  2*Gz*Pz*C + 2*Gz*P*Cz + 2*G*Pz*Cz
146
147    if return_components:
148        return grad_2_G_unmixed_xx, grad_2_G_unmixed_yy, grad_2_G_unmixed_zz, G,P,C,Ga,Pa,Ca
149    return grad_2_G_unmixed_xx, grad_2_G_unmixed_yy, grad_2_G_unmixed_zz
150
151
152def get_G_second_mixed(points:Tensor, scatterer:Mesh, board:Tensor|None=None, return_components:bool=False, k=Constants.k) -> tuple[Tensor, Tensor, Tensor]:
153
154    #Bk.3 Pg.81
155    
156    areas = get_areas(scatterer).unsqueeze(1)
157    centres = get_centres_as_points(scatterer)
158    normals = get_normals_as_points(scatterer).unsqueeze(2)
159
160
161    grad_Gx, grad_Gy, grad_Gz, G,P,C,Ga,Pa,Ca = get_G_partial(points, scatterer, board, True, k=k)
162
163    points  = points.unsqueeze(3)  # [B, 3, N, 1]
164    centres = centres.unsqueeze(2)  # [B, 3, 1, M]
165
166    diff = points - centres    
167    diff_square = diff**2
168    distances = torch.sqrt(torch.sum(diff_square, 1))
169    distances_square = distances ** 2
170    distances_cube = distances ** 3
171    
172    kd = k * distances
173    ikd = 1j * kd
174
175    dx = diff[:,0,:,:]
176    dy = diff[:,1,:,:]
177    dz = diff[:,2,:,:]
178
179    # print(dx.shape, dy.shape, dz.shape, distances_cube.shape)
180    dxy = -dx*dy / distances_cube
181    dxz = -dx*dz / distances_cube
182    dyz = -dy*dz / distances_cube
183
184    # print(dx.shape, distances.shape)
185    dax = dx / distances
186    day = dy / distances
187    daz = dz / distances
188
189    exp_ikd = torch.exp(ikd)
190
191
192
193    # print(areas.shape, exp_ikd.shape, day.shape, dax.shape, distances.shape, distances_cube.shape)
194    Gxy = (-areas * exp_ikd * (day * dax * (kd**2 + 2*ikd - 2) + distances * dxy * (1 - ikd))) / (4*Constants.pi * distances_cube)
195    Gxz = (-areas * exp_ikd * (daz * dax * (kd**2 + 2*ikd - 2) + distances * dxz * (1 - ikd))) / (4*Constants.pi * distances_cube)
196    Gyz = (-areas * exp_ikd * (day * daz * (kd**2 + 2*ikd - 2) + distances * dyz * (1 - ikd))) / (4*Constants.pi * distances_cube)
197
198    nx = normals[:,0,:]
199    ny = normals[:,1,:]
200    nz = normals[:,2,:]
201
202    # print(nx.shape, dx.shape, ny.shape, dy.shape, nz.shape, dz.shape)
203    nd = nx * dx + ny * dy + nz * dz
204
205    Cxy = (-nx*dy - ny*dx + (3*nd*dx*dy/distances_cube)) / distances_square
206    Cxz = (-nx*dz - nz*dx + (3*nd*dx*dz/distances_cube)) / distances_square
207    Cyz = (-nz*dy - ny*dz + (3*nd*dz*dy/distances_cube)) / distances_square
208
209
210    Pxy = (distances * dxy - 2 * day * dax) / distances_cube
211    Pxz = (distances * dxz - 2 * daz * dax) / distances_cube
212    Pyz = (distances * dyz - 2 * day * daz) / distances_cube
213
214    
215    G = G[:,0,:]
216    P = P[:,0,:] 
217    # C = C.unsqueeze(1).expand(-1,3,-1,-1)
218
219    Gx = Ga[:,0,:]
220    Gy = Ga[:,1,:]
221    Gz = Ga[:,2,:]
222
223    Px = Pa[:,0,:]
224    Py = Pa[:,1,:]
225    Pz = Pa[:,2,:]
226
227    Cx = Ca[:,0,:]
228    Cy = Ca[:,1,:]
229    Cz = Ca[:,2,:]
230
231    # grad_2_G_mixed_xy = Gxy*P*C + Gy*Px*C + Gy*P*Cx + Gx*Py*C + G*Pxy*C + G*Py*Cx + Gx*P*Cy + G*Px*Cy + G*P*Cxy
232    # grad_2_G_mixed_xz = Gxz*P*C + Gz*Px*C + Gz*P*Cx + Gx*Pz*C + G*Pxz*C + G*Pz*Cx + Gx*P*Cz + G*Px*Cz + G*P*Cxz
233    # grad_2_G_mixed_yz = Gyz*P*C + Gy*Pz*C + Gy*P*Cz + Gz*Py*C + G*Pyz*C + G*Py*Cz + Gz*P*Cy + G*Pz*Cy + G*P*Cyz
234
235
236    grad_2_G_mixed_xy = Gxy*P*C + Gx*Py*C + G*Pxy*C + G*Py*Cx + Gx*P*Cy + G*P*Cxy
237    grad_2_G_mixed_xz = Gxz*P*C + Gx*Pz*C + G*Pxz*C + G*Pz*Cx + Gx*P*Cz + G*P*Cxz
238    grad_2_G_mixed_yz = Gyz*P*C + Gz*Py*C + G*Pyz*C + G*Py*Cz + Gz*P*Cy + G*P*Cyz
239
240
241
242    return grad_2_G_mixed_xy, grad_2_G_mixed_xz, grad_2_G_mixed_yz
243
244
245def BEM_forward_model_second_derivative_unmixed(points:Tensor, scatterer:Mesh, transducers:Tensor=None, use_cache_H:bool=True, k=Constants.k,
246                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False,
247                           path:str="Media", p_ref=Constants.P_ref,internal_points = None, transducer_radius=Constants.radius):
248                           
249    
250    
251    if transducers is None:
252        transducers = TRANSDUCERS
253
254    if H is None:
255        H = get_cache_or_compute_H(scatterer,transducers,use_cache_H, path, print_lines,p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius)
256
257    Fxx, Fyy, Fzz = forward_model_second_derivative_unmixed(points,transducers=transducers,p_ref=p_ref, k=k, transducer_radius=transducer_radius)
258    Gxx, Gyy, Gzz = get_G_second_unmixed(points, scatterer, transducers, k=k)
259
260    Exx = Fxx + Gxx@H
261    Eyy = Fyy + Gyy@H
262    Ezz = Fzz + Gzz@H
263
264    return Exx, Eyy, Ezz
265
266def BEM_forward_model_second_derivative_mixed(points:Tensor, scatterer:Mesh, transducers:Tensor|Mesh=None, use_cache_H:bool=True, k=Constants.k,
267                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False, 
268                           path:str="Media", p_ref=Constants.P_ref,internal_points = None, transducer_radius=Constants.radius):
269    
270       
271    if transducers is None:
272        transducers = TRANSDUCERS
273
274    if H is None:
275        H = get_cache_or_compute_H(scatterer,transducers,use_cache_H, path, print_lines,p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius)
276
277    Fxy, Fxz, Fyz = forward_model_second_derivative_mixed(points,transducers=transducers, k=k, transducer_radius=transducer_radius, p_ref=p_ref)
278    Gxy, Gxz, Gyz = get_G_second_mixed(points, scatterer, transducers, k=k)
279
280
281    Exy = Fxy + Gxy@H
282    Exz = Fxz + Gxz@H
283    Eyz = Fyz + Gyz@H
284
285    return Exy, Exz, Eyz
286
287def BEM_laplacian(points:Tensor, scatterer:Mesh, transducers:Tensor|Mesh=None, use_cache_H:bool=True, k=Constants.k,
288                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False, 
289                           path:str="Media", p_ref=Constants.P_ref,internal_points = None, transducer_radius=Constants.radius):
290    
291    '''
292    Computes the laplacian of pressure at points given a hologram
293
294    :param activations: Transducer hologram
295    :param points: Points to propagate to
296    :param scatterer: The mesh used (as a `vedo` `mesh` object)
297    :param board: Transducers to use, if `None` then uses `acoustools.Utilities.TOP_BOARD` 
298    :param H: Precomputed H - if None H will be computed
299    :param E: Precomputed E - if None E will be computed
300    :param path: path to folder containing `BEMCache/ `
301    :param use_cache_H: If True uses the cache system to load and save the H matrix. Default `True`
302    :param print_lines: if true prints messages detaling progress
303    :param k: wavenumber
304    :param internal_points: The internal points to use for CHIEF based BEM
305
306    :return pressure: laplacian at points
307    '''
308    
309    Exx, Eyy, Ezz = BEM_forward_model_second_derivative_unmixed(points=points, scatterer=scatterer, transducers=transducers, use_cache_H=use_cache_H, k=k, print_lines=print_lines, H=H, return_components=return_components, path=path, p_ref=p_ref, internal_points=internal_points, transducer_radius=transducer_radius)
310
311    return Exx + Eyy + Ezz
312
313
314
315def get_G_partial(points:Tensor, scatterer:Mesh, board:Tensor|None=None, return_components:bool=False, k=Constants.k) -> tuple[Tensor, Tensor, Tensor]:
316    '''
317    Computes gradient of the G matrix in BEM \n
318    :param points: Points to propagate to
319    :param scatterer: The mesh used (as a `vedo` `mesh` object)
320    :param board: Ignored
321    :param return_components: if true will return the subparts used to compute
322    :return: Gradient of the G matrix in BEM
323    '''
324    #Bk3. Pg. 26
325    # if board is None:
326    #     board = TRANSDUCERS
327
328    areas = get_areas(scatterer)
329    centres = get_centres_as_points(scatterer)
330    normals = get_normals_as_points(scatterer)
331
332
333    N = points.shape[2]
334    M = centres.shape[2]
335
336
337    # points = points.unsqueeze(3).expand(-1,-1,-1,M)
338    # centres = centres.unsqueeze(2).expand(-1,-1,N,-1)
339    points  = points.unsqueeze(3)  # [B, 3, N, 1]
340    centres = centres.unsqueeze(2)  # [B, 3, 1, M]
341
342    diff = (points - centres)
343    diff_square = diff**2
344    distances = torch.sqrt(torch.sum(diff_square, 1))
345    distances_expanded = distances.unsqueeze(1)#.expand((1,3,N,M))
346    distances_expanded_square = distances_expanded**2
347    distances_expanded_cube = distances_expanded**3
348
349    # G  =  e^(ikd) / 4pi d
350    G = areas * torch.exp(1j * k * distances_expanded) / (4*3.1415*distances_expanded)
351
352    #Ga =  [i*da * e^{ikd} * (kd+i) / 4pi d^2]
353
354    #d = distance
355    #da = -(at - a)^2 / d
356    da = diff / distances_expanded
357    kd = k * distances_expanded
358    phase = torch.exp(1j*kd)
359    Ga =  areas * ( (1j*da*phase * (kd + 1j))/ (4*3.1415*distances_expanded_square))
360
361    #P = (ik - 1/d)
362    P = (1j*k - 1/distances_expanded)
363    #Pa = da / d^2 = (diff / d^2) /d
364    Pa = diff / distances_expanded_cube
365
366    #C = (diff \cdot normals) / distances
367
368    nx = normals[:,0]
369    ny = normals[:,1]
370    nz = normals[:,2]
371
372    dx = diff[:,0,:]
373    dy = diff[:,1,:]
374    dz = diff[:,2,:]
375
376    n_dot_d = nx*dx + ny*dy + nz*dz
377
378    C = (n_dot_d) / distances
379
380
381    distance_square = distances**2
382
383
384    Cx = 1/distance_square * (nx * distances - (n_dot_d * dx) / distances)
385    Cy = 1/distance_square * (ny * distances - (n_dot_d * dy) / distances)
386    Cz = 1/distance_square * (nz * distances - (n_dot_d * dz) / distances)
387
388    Cx.unsqueeze_(1)
389    Cy.unsqueeze_(1)
390    Cz.unsqueeze_(1)
391
392    Ca = torch.cat([Cx, Cy, Cz],axis=1)
393
394    grad_G = Ga*P*C + G*P*Ca + G*Pa*C
395
396    grad_G =  -grad_G.to(DTYPE)
397
398    if return_components:
399        return grad_G[:,0,:], grad_G[:,1,:], grad_G[:,2,:], G,P,C,Ga,Pa, Ca
400    
401    return grad_G[:,0,:], grad_G[:,1,:], grad_G[:,2,:]

def BEM_forward_model_grad( points: torch.Tensor, scatterer: vedo.mesh.Mesh, transducers: torch.Tensor = None, use_cache_H: bool = True, print_lines: bool = False, H: torch.Tensor | None = None, return_components: bool = False, k=732.7329804081634, path: str = 'Media', p_ref=3.4000000000000004, internal_points=None, transducer_radius=0.0045) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor] | tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]: View Source

15def BEM_forward_model_grad(points:Tensor, scatterer:Mesh, transducers:Tensor=None, use_cache_H:bool=True, 
16                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False,k=Constants.k,
17                           path:str="Media", p_ref=Constants.P_ref, internal_points = None, transducer_radius=Constants.radius) -> tuple[Tensor, Tensor, Tensor] | tuple[Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor]:
18    '''
19    Computes the gradient of the forward propagation for BEM\n
20    :param scatterer: The mesh used (as a `vedo` `mesh` object)
21    :param transducers: Transducers to use, if `None` uses `acoustools.Utilities.TRANSDUCERS`
22    :param use_cache_H_grad: If true uses the cache system, otherwise computes `H` and does not save it
23    :param print_lines: if true prints messages detaling progress
24    :param H: Precomputed `H` - if `None` `H` will be computed
25    :param return_components: if true will return the subparts used to compute
26    :param path: path to folder containing `BEMCache/` 
27    :return: Ex, Ey, Ez
28    '''
29    if transducers is None:
30        transducers = TRANSDUCERS
31
32    B = points.shape[0]
33    if H is None:
34        H = get_cache_or_compute_H(scatterer,transducers,use_cache_H, path, print_lines,p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius)
35    
36    Fx, Fy, Fz  = forward_model_grad(points, transducers,p_ref=p_ref, k=k, transducer_radius=transducer_radius)
37    Gx, Gy, Gz = get_G_partial(points, scatterer, transducers, k=k)
38
39
40    Ex = Fx + Gx@H 
41    Ey = Fy + Gy@H
42    Ez = Fz + Gz@H
43
44
45    if return_components:
46        return Ex.to(DTYPE), Ey.to(DTYPE), Ez.to(DTYPE), Fx, Fy, Fz, Gx, Gy, Gz, H
47    else:
48        return Ex.to(DTYPE), Ey.to(DTYPE), Ez.to(DTYPE)

Computes the gradient of the forward propagation for BEM

Parameters

scatterer: The mesh used (as a vedo mesh object)
transducers: Transducers to use, if None uses acoustools.Utilities.TRANSDUCERS
use_cache_H_grad: If true uses the cache system, otherwise computes H and does not save it
print_lines: if true prints messages detaling progress
H: Precomputed H - if None H will be computed
return_components: if true will return the subparts used to compute
path: path to folder containing BEMCache/

Returns

Ex, Ey, Ez

def get_G_second_unmixed( points: torch.Tensor, scatterer: vedo.mesh.Mesh, board: torch.Tensor | None = None, return_components: bool = False, k=732.7329804081634) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]: View Source

 51def get_G_second_unmixed(points:Tensor, scatterer:Mesh, board:Tensor|None=None, return_components:bool=False, k=Constants.k) -> tuple[Tensor, Tensor, Tensor]:
 52
 53    #Bk.3 Pg.80
 54    #P is much higher than the others?
 55
 56    areas = get_areas(scatterer)
 57    centres = get_centres_as_points(scatterer)
 58    normals = get_normals_as_points(scatterer).unsqueeze(2)
 59
 60    N = points.shape[2]
 61    M = centres.shape[2]
 62
 63    points.requires_grad_()
 64    grad_Gx, grad_Gy, grad_Gz, G,P,C,Ga,Pa,Ca = get_G_partial(points, scatterer, board, True, k=k)
 65    # exit()
 66
 67    points  = points.unsqueeze(3)  # [B, 3, N, 1]
 68    centres = centres.unsqueeze(2)  # [B, 3, 1, M]
 69    
 70
 71    diff = points - centres    
 72    diff_square = diff**2
 73    diff_square_sum = torch.sum(diff_square, 1)
 74    distances = torch.sqrt(diff_square_sum)
 75    distances_cube = distances ** 3
 76    distances_five = distances ** 5
 77
 78    distances_expanded = distances.unsqueeze(1).expand((1,3,N,M))
 79    distances_expanded_square = distances_expanded**2
 80    distances_expanded_cube = distances_expanded**3
 81    distances_expanded_four = distances_expanded**4
 82    distances_expanded_five = distances_expanded**5
 83
 84    kd = k * distances_expanded
 85
 86    exp_ikd = torch.exp(1j*kd)
 87    
 88    Gaa = areas/(4*Constants.pi * 1j) *exp_ikd * ((1j*k**2 * diff_square)/ distances_expanded_cube + (1j*k)/distances_expanded_square - 1/distances_expanded_cube - (3*1j*k*diff_square)/distances_expanded_four + (3*diff_square)/distances_expanded_five)
 89    Gxx = Gaa[:,0]
 90    Gyy = Gaa[:,1]
 91    Gzz = Gaa[:,2]
 92    
 93    #P = ik - 1/d
 94    # Paa = (distances_expanded * daa - 2*da**2) / distances_expanded_cube
 95    # Paa = (3 * diff_square - distances_expanded_square) / distances_expanded_five
 96    Paa = 1/distances_expanded_cube - (3*diff_square) / distances_expanded_five  #Equal to lines above
 97    Pxx = Paa[:,0]
 98    Pyy = Paa[:,1]
 99    Pzz = Paa[:,2]
100
101    #C = (diff dot c )/ d
102
103    dx = diff[:,0,:,:]
104    dy = diff[:,1,:,:]
105    dz = diff[:,2,:,:]
106
107
108    nx = normals[:,0,:]
109    ny = normals[:,1,:]
110    nz = normals[:,2,:]
111
112    nd = nx * dx + ny * dy + nz * dz
113
114    Cxx = (2 * nx * dx)/distances_cube - nd/distances_cube + (3*nd*dx**2)/distances_five
115    Cyy = (2 * ny * dy)/distances_cube - nd/distances_cube + (3*nd*dy**2)/distances_five
116    Czz = (2 * nz * dz)/distances_cube - nd/distances_cube + (3*nd*dz**2)/distances_five
117
118    # Caa = torch.stack([Cxx, Cyy, Czz], dim=1)
119
120    # Cxx = (((3 * dax**2)/distances_five) - (1/distances_cube)) * nd - (2*nx*dx) / distances_cube
121    # Cyy = (((3 * day**2)/distances_five) - (1/distances_cube)) * nd - (2*ny*dy) / distances_cube 
122    # Czz = (((3 * daz**2)/distances_five) - (1/distances_cube)) * nd - (2*nz*dz) / distances_cube 
123    # Caa = torch.stack([Cxx, Cyy, Czz], dim=1) 
124
125    # grad_2_G_unmixed_old = 2 * Ca * (Ga * P + G * Pa) + C * (2 * Ga*Pa + Gaa * P + G * Paa) + G*P*Caa
126    # grad_2_G_unmixed = C * (2*Ga * Pa + Gaa * P + G*Paa) + P * (2*Ga * Ca + G*Caa) + 2*G * Pa * Ca
127
128    Cx = Ca[:,0]
129    Cy = Ca[:,1]
130    Cz = Ca[:,2]
131
132    Px = Pa[:,0]
133    Py = Pa[:,1]
134    Pz = Pa[:,2]
135
136    Gx = Ga[:,0]
137    Gy = Ga[:,1]
138    Gz = Ga[:,2]
139
140    G = G[:,0,:]
141    P = P[:,0,:] 
142    
143
144
145    grad_2_G_unmixed_xx = Gxx * P * C + G*Pxx*C + G*P*Cxx+  2*Gx*Px*C + 2*Gx*P*Cx + 2*G*Px*Cx
146    grad_2_G_unmixed_yy = Gyy * P * C + G*Pyy*C + G*P*Cyy+  2*Gy*Py*C + 2*Gy*P*Cy + 2*G*Py*Cy
147    grad_2_G_unmixed_zz = Gzz * P * C + G*Pzz*C + G*P*Czz+  2*Gz*Pz*C + 2*Gz*P*Cz + 2*G*Pz*Cz
148
149    if return_components:
150        return grad_2_G_unmixed_xx, grad_2_G_unmixed_yy, grad_2_G_unmixed_zz, G,P,C,Ga,Pa,Ca
151    return grad_2_G_unmixed_xx, grad_2_G_unmixed_yy, grad_2_G_unmixed_zz

def get_G_second_mixed( points: torch.Tensor, scatterer: vedo.mesh.Mesh, board: torch.Tensor | None = None, return_components: bool = False, k=732.7329804081634) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]: View Source

154def get_G_second_mixed(points:Tensor, scatterer:Mesh, board:Tensor|None=None, return_components:bool=False, k=Constants.k) -> tuple[Tensor, Tensor, Tensor]:
155
156    #Bk.3 Pg.81
157    
158    areas = get_areas(scatterer).unsqueeze(1)
159    centres = get_centres_as_points(scatterer)
160    normals = get_normals_as_points(scatterer).unsqueeze(2)
161
162
163    grad_Gx, grad_Gy, grad_Gz, G,P,C,Ga,Pa,Ca = get_G_partial(points, scatterer, board, True, k=k)
164
165    points  = points.unsqueeze(3)  # [B, 3, N, 1]
166    centres = centres.unsqueeze(2)  # [B, 3, 1, M]
167
168    diff = points - centres    
169    diff_square = diff**2
170    distances = torch.sqrt(torch.sum(diff_square, 1))
171    distances_square = distances ** 2
172    distances_cube = distances ** 3
173    
174    kd = k * distances
175    ikd = 1j * kd
176
177    dx = diff[:,0,:,:]
178    dy = diff[:,1,:,:]
179    dz = diff[:,2,:,:]
180
181    # print(dx.shape, dy.shape, dz.shape, distances_cube.shape)
182    dxy = -dx*dy / distances_cube
183    dxz = -dx*dz / distances_cube
184    dyz = -dy*dz / distances_cube
185
186    # print(dx.shape, distances.shape)
187    dax = dx / distances
188    day = dy / distances
189    daz = dz / distances
190
191    exp_ikd = torch.exp(ikd)
192
193
194
195    # print(areas.shape, exp_ikd.shape, day.shape, dax.shape, distances.shape, distances_cube.shape)
196    Gxy = (-areas * exp_ikd * (day * dax * (kd**2 + 2*ikd - 2) + distances * dxy * (1 - ikd))) / (4*Constants.pi * distances_cube)
197    Gxz = (-areas * exp_ikd * (daz * dax * (kd**2 + 2*ikd - 2) + distances * dxz * (1 - ikd))) / (4*Constants.pi * distances_cube)
198    Gyz = (-areas * exp_ikd * (day * daz * (kd**2 + 2*ikd - 2) + distances * dyz * (1 - ikd))) / (4*Constants.pi * distances_cube)
199
200    nx = normals[:,0,:]
201    ny = normals[:,1,:]
202    nz = normals[:,2,:]
203
204    # print(nx.shape, dx.shape, ny.shape, dy.shape, nz.shape, dz.shape)
205    nd = nx * dx + ny * dy + nz * dz
206
207    Cxy = (-nx*dy - ny*dx + (3*nd*dx*dy/distances_cube)) / distances_square
208    Cxz = (-nx*dz - nz*dx + (3*nd*dx*dz/distances_cube)) / distances_square
209    Cyz = (-nz*dy - ny*dz + (3*nd*dz*dy/distances_cube)) / distances_square
210
211
212    Pxy = (distances * dxy - 2 * day * dax) / distances_cube
213    Pxz = (distances * dxz - 2 * daz * dax) / distances_cube
214    Pyz = (distances * dyz - 2 * day * daz) / distances_cube
215
216    
217    G = G[:,0,:]
218    P = P[:,0,:] 
219    # C = C.unsqueeze(1).expand(-1,3,-1,-1)
220
221    Gx = Ga[:,0,:]
222    Gy = Ga[:,1,:]
223    Gz = Ga[:,2,:]
224
225    Px = Pa[:,0,:]
226    Py = Pa[:,1,:]
227    Pz = Pa[:,2,:]
228
229    Cx = Ca[:,0,:]
230    Cy = Ca[:,1,:]
231    Cz = Ca[:,2,:]
232
233    # grad_2_G_mixed_xy = Gxy*P*C + Gy*Px*C + Gy*P*Cx + Gx*Py*C + G*Pxy*C + G*Py*Cx + Gx*P*Cy + G*Px*Cy + G*P*Cxy
234    # grad_2_G_mixed_xz = Gxz*P*C + Gz*Px*C + Gz*P*Cx + Gx*Pz*C + G*Pxz*C + G*Pz*Cx + Gx*P*Cz + G*Px*Cz + G*P*Cxz
235    # grad_2_G_mixed_yz = Gyz*P*C + Gy*Pz*C + Gy*P*Cz + Gz*Py*C + G*Pyz*C + G*Py*Cz + Gz*P*Cy + G*Pz*Cy + G*P*Cyz
236
237
238    grad_2_G_mixed_xy = Gxy*P*C + Gx*Py*C + G*Pxy*C + G*Py*Cx + Gx*P*Cy + G*P*Cxy
239    grad_2_G_mixed_xz = Gxz*P*C + Gx*Pz*C + G*Pxz*C + G*Pz*Cx + Gx*P*Cz + G*P*Cxz
240    grad_2_G_mixed_yz = Gyz*P*C + Gz*Py*C + G*Pyz*C + G*Py*Cz + Gz*P*Cy + G*P*Cyz
241
242
243
244    return grad_2_G_mixed_xy, grad_2_G_mixed_xz, grad_2_G_mixed_yz

def BEM_forward_model_second_derivative_unmixed( points: torch.Tensor, scatterer: vedo.mesh.Mesh, transducers: torch.Tensor = None, use_cache_H: bool = True, k=732.7329804081634, print_lines: bool = False, H: torch.Tensor | None = None, return_components: bool = False, path: str = 'Media', p_ref=3.4000000000000004, internal_points=None, transducer_radius=0.0045): View Source

247def BEM_forward_model_second_derivative_unmixed(points:Tensor, scatterer:Mesh, transducers:Tensor=None, use_cache_H:bool=True, k=Constants.k,
248                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False,
249                           path:str="Media", p_ref=Constants.P_ref,internal_points = None, transducer_radius=Constants.radius):
250                           
251    
252    
253    if transducers is None:
254        transducers = TRANSDUCERS
255
256    if H is None:
257        H = get_cache_or_compute_H(scatterer,transducers,use_cache_H, path, print_lines,p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius)
258
259    Fxx, Fyy, Fzz = forward_model_second_derivative_unmixed(points,transducers=transducers,p_ref=p_ref, k=k, transducer_radius=transducer_radius)
260    Gxx, Gyy, Gzz = get_G_second_unmixed(points, scatterer, transducers, k=k)
261
262    Exx = Fxx + Gxx@H
263    Eyy = Fyy + Gyy@H
264    Ezz = Fzz + Gzz@H
265
266    return Exx, Eyy, Ezz

def BEM_forward_model_second_derivative_mixed( points: torch.Tensor, scatterer: vedo.mesh.Mesh, transducers: torch.Tensor | vedo.mesh.Mesh = None, use_cache_H: bool = True, k=732.7329804081634, print_lines: bool = False, H: torch.Tensor | None = None, return_components: bool = False, path: str = 'Media', p_ref=3.4000000000000004, internal_points=None, transducer_radius=0.0045): View Source

268def BEM_forward_model_second_derivative_mixed(points:Tensor, scatterer:Mesh, transducers:Tensor|Mesh=None, use_cache_H:bool=True, k=Constants.k,
269                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False, 
270                           path:str="Media", p_ref=Constants.P_ref,internal_points = None, transducer_radius=Constants.radius):
271    
272       
273    if transducers is None:
274        transducers = TRANSDUCERS
275
276    if H is None:
277        H = get_cache_or_compute_H(scatterer,transducers,use_cache_H, path, print_lines,p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius)
278
279    Fxy, Fxz, Fyz = forward_model_second_derivative_mixed(points,transducers=transducers, k=k, transducer_radius=transducer_radius, p_ref=p_ref)
280    Gxy, Gxz, Gyz = get_G_second_mixed(points, scatterer, transducers, k=k)
281
282
283    Exy = Fxy + Gxy@H
284    Exz = Fxz + Gxz@H
285    Eyz = Fyz + Gyz@H
286
287    return Exy, Exz, Eyz

def BEM_laplacian( points: torch.Tensor, scatterer: vedo.mesh.Mesh, transducers: torch.Tensor | vedo.mesh.Mesh = None, use_cache_H: bool = True, k=732.7329804081634, print_lines: bool = False, H: torch.Tensor | None = None, return_components: bool = False, path: str = 'Media', p_ref=3.4000000000000004, internal_points=None, transducer_radius=0.0045): View Source

289def BEM_laplacian(points:Tensor, scatterer:Mesh, transducers:Tensor|Mesh=None, use_cache_H:bool=True, k=Constants.k,
290                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False, 
291                           path:str="Media", p_ref=Constants.P_ref,internal_points = None, transducer_radius=Constants.radius):
292    
293    '''
294    Computes the laplacian of pressure at points given a hologram
295
296    :param activations: Transducer hologram
297    :param points: Points to propagate to
298    :param scatterer: The mesh used (as a `vedo` `mesh` object)
299    :param board: Transducers to use, if `None` then uses `acoustools.Utilities.TOP_BOARD` 
300    :param H: Precomputed H - if None H will be computed
301    :param E: Precomputed E - if None E will be computed
302    :param path: path to folder containing `BEMCache/ `
303    :param use_cache_H: If True uses the cache system to load and save the H matrix. Default `True`
304    :param print_lines: if true prints messages detaling progress
305    :param k: wavenumber
306    :param internal_points: The internal points to use for CHIEF based BEM
307
308    :return pressure: laplacian at points
309    '''
310    
311    Exx, Eyy, Ezz = BEM_forward_model_second_derivative_unmixed(points=points, scatterer=scatterer, transducers=transducers, use_cache_H=use_cache_H, k=k, print_lines=print_lines, H=H, return_components=return_components, path=path, p_ref=p_ref, internal_points=internal_points, transducer_radius=transducer_radius)
312
313    return Exx + Eyy + Ezz

Computes the laplacian of pressure at points given a hologram

Parameters

activations: Transducer hologram
points: Points to propagate to
scatterer: The mesh used (as a vedo mesh object)
board: Transducers to use, if None then uses acoustools.Utilities.TOP_BOARD
H: Precomputed H - if None H will be computed
E: Precomputed E - if None E will be computed
path: path to folder containing BEMCache/
use_cache_H: If True uses the cache system to load and save the H matrix. Default True
print_lines: if true prints messages detaling progress
k: wavenumber
internal_points: The internal points to use for CHIEF based BEM

Returns

laplacian at points

def get_G_partial( points: torch.Tensor, scatterer: vedo.mesh.Mesh, board: torch.Tensor | None = None, return_components: bool = False, k=732.7329804081634) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]: View Source

317def get_G_partial(points:Tensor, scatterer:Mesh, board:Tensor|None=None, return_components:bool=False, k=Constants.k) -> tuple[Tensor, Tensor, Tensor]:
318    '''
319    Computes gradient of the G matrix in BEM \n
320    :param points: Points to propagate to
321    :param scatterer: The mesh used (as a `vedo` `mesh` object)
322    :param board: Ignored
323    :param return_components: if true will return the subparts used to compute
324    :return: Gradient of the G matrix in BEM
325    '''
326    #Bk3. Pg. 26
327    # if board is None:
328    #     board = TRANSDUCERS
329
330    areas = get_areas(scatterer)
331    centres = get_centres_as_points(scatterer)
332    normals = get_normals_as_points(scatterer)
333
334
335    N = points.shape[2]
336    M = centres.shape[2]
337
338
339    # points = points.unsqueeze(3).expand(-1,-1,-1,M)
340    # centres = centres.unsqueeze(2).expand(-1,-1,N,-1)
341    points  = points.unsqueeze(3)  # [B, 3, N, 1]
342    centres = centres.unsqueeze(2)  # [B, 3, 1, M]
343
344    diff = (points - centres)
345    diff_square = diff**2
346    distances = torch.sqrt(torch.sum(diff_square, 1))
347    distances_expanded = distances.unsqueeze(1)#.expand((1,3,N,M))
348    distances_expanded_square = distances_expanded**2
349    distances_expanded_cube = distances_expanded**3
350
351    # G  =  e^(ikd) / 4pi d
352    G = areas * torch.exp(1j * k * distances_expanded) / (4*3.1415*distances_expanded)
353
354    #Ga =  [i*da * e^{ikd} * (kd+i) / 4pi d^2]
355
356    #d = distance
357    #da = -(at - a)^2 / d
358    da = diff / distances_expanded
359    kd = k * distances_expanded
360    phase = torch.exp(1j*kd)
361    Ga =  areas * ( (1j*da*phase * (kd + 1j))/ (4*3.1415*distances_expanded_square))
362
363    #P = (ik - 1/d)
364    P = (1j*k - 1/distances_expanded)
365    #Pa = da / d^2 = (diff / d^2) /d
366    Pa = diff / distances_expanded_cube
367
368    #C = (diff \cdot normals) / distances
369
370    nx = normals[:,0]
371    ny = normals[:,1]
372    nz = normals[:,2]
373
374    dx = diff[:,0,:]
375    dy = diff[:,1,:]
376    dz = diff[:,2,:]
377
378    n_dot_d = nx*dx + ny*dy + nz*dz
379
380    C = (n_dot_d) / distances
381
382
383    distance_square = distances**2
384
385
386    Cx = 1/distance_square * (nx * distances - (n_dot_d * dx) / distances)
387    Cy = 1/distance_square * (ny * distances - (n_dot_d * dy) / distances)
388    Cz = 1/distance_square * (nz * distances - (n_dot_d * dz) / distances)
389
390    Cx.unsqueeze_(1)
391    Cy.unsqueeze_(1)
392    Cz.unsqueeze_(1)
393
394    Ca = torch.cat([Cx, Cy, Cz],axis=1)
395
396    grad_G = Ga*P*C + G*P*Ca + G*Pa*C
397
398    grad_G =  -grad_G.to(DTYPE)
399
400    if return_components:
401        return grad_G[:,0,:], grad_G[:,1,:], grad_G[:,2,:], G,P,C,Ga,Pa, Ca
402    
403    return grad_G[:,0,:], grad_G[:,1,:], grad_G[:,2,:]

Computes gradient of the G matrix in BEM

Parameters

points: Points to propagate to
scatterer: The mesh used (as a vedo mesh object)
board: Ignored
return_components: if true will return the subparts used to compute

Returns

Gradient of the G matrix in BEM