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

def BEM_forward_model_grad( points: torch.Tensor, scatterer: vedo.mesh.core.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, transducer_norms=None) -> 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, transducer_norms=None) -> 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, norms=transducer_norms)
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.core.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    ik = 1j * k
 85    ikd = ik * distances_expanded
 86
 87    exp_ikd = torch.exp(ikd)
 88    
 89    Gaa = areas/(4*Constants.pi * 1j) *exp_ikd * ((1j*k**2 * diff_square)/ distances_expanded_cube + (ik)/distances_expanded_square - 1/distances_expanded_cube - (3*ik*diff_square)/distances_expanded_four + (3*diff_square)/distances_expanded_five)
 90    Gxx = Gaa[:,0]
 91    Gyy = Gaa[:,1]
 92    Gzz = Gaa[:,2]
 93    
 94    #P = ik - 1/d
 95    # Paa = (distances_expanded * daa - 2*da**2) / distances_expanded_cube
 96    # Paa = (3 * diff_square - distances_expanded_square) / distances_expanded_five
 97    Paa = 1/distances_expanded_cube - (3*diff_square) / distances_expanded_five  #Equal to lines above
 98    Pxx = Paa[:,0]
 99    Pyy = Paa[:,1]
100    Pzz = Paa[:,2]
101
102    #C = (diff dot c )/ d
103
104    dx = diff[:,0,:,:]
105    dy = diff[:,1,:,:]
106    dz = diff[:,2,:,:]
107
108
109    nx = normals[:,0,:]
110    ny = normals[:,1,:]
111    nz = normals[:,2,:]
112
113    nd = nx * dx + ny * dy + nz * dz
114
115    Cxx = (2 * nx * dx)/distances_cube - nd/distances_cube + (3*nd*dx**2)/distances_five
116    Cyy = (2 * ny * dy)/distances_cube - nd/distances_cube + (3*nd*dy**2)/distances_five
117    Czz = (2 * nz * dz)/distances_cube - nd/distances_cube + (3*nd*dz**2)/distances_five
118
119    # Caa = torch.stack([Cxx, Cyy, Czz], dim=1)
120
121    # Cxx = (((3 * dax**2)/distances_five) - (1/distances_cube)) * nd - (2*nx*dx) / distances_cube
122    # Cyy = (((3 * day**2)/distances_five) - (1/distances_cube)) * nd - (2*ny*dy) / distances_cube 
123    # Czz = (((3 * daz**2)/distances_five) - (1/distances_cube)) * nd - (2*nz*dz) / distances_cube 
124    # Caa = torch.stack([Cxx, Cyy, Czz], dim=1) 
125
126    # grad_2_G_unmixed_old = 2 * Ca * (Ga * P + G * Pa) + C * (2 * Ga*Pa + Gaa * P + G * Paa) + G*P*Caa
127    # grad_2_G_unmixed = C * (2*Ga * Pa + Gaa * P + G*Paa) + P * (2*Ga * Ca + G*Caa) + 2*G * Pa * Ca
128
129    Cx = Ca[:,0]
130    Cy = Ca[:,1]
131    Cz = Ca[:,2]
132
133    Px = Pa[:,0]
134    Py = Pa[:,1]
135    Pz = Pa[:,2]
136
137    Gx = Ga[:,0]
138    Gy = Ga[:,1]
139    Gz = Ga[:,2]
140
141    G = G[:,0,:]
142    P = P[:,0,:] 
143    
144
145
146    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
147    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
148    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
149
150    if return_components:
151        return grad_2_G_unmixed_xx, grad_2_G_unmixed_yy, grad_2_G_unmixed_zz, G,P,C,Ga,Pa,Ca
152    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.core.Mesh, board: torch.Tensor | None = None, return_components: bool = False, k=732.7329804081634) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]: View Source

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

251def BEM_forward_model_second_derivative_unmixed(points:Tensor, scatterer:Mesh, transducers:Tensor=None, use_cache_H:bool=True, k=Constants.k,
252                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False,
253                           path:str="Media", p_ref=Constants.P_ref,internal_points = None, transducer_radius=Constants.radius, transducer_norms=None):
254                           
255    
256    
257    if transducers is None:
258        transducers = TRANSDUCERS
259
260    if H is None:
261        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, norms=transducer_norms)
262
263    Fxx, Fyy, Fzz = forward_model_second_derivative_unmixed(points,transducers=transducers,p_ref=p_ref, k=k, transducer_radius=transducer_radius, transducer_norms=transducer_norms)
264    Gxx, Gyy, Gzz = get_G_second_unmixed(points, scatterer, transducers, k=k)
265
266    Exx = Fxx + Gxx@H
267    Eyy = Fyy + Gyy@H
268    Ezz = Fzz + Gzz@H
269
270    return Exx, Eyy, Ezz

def BEM_forward_model_second_derivative_mixed( points: torch.Tensor, scatterer: vedo.mesh.core.Mesh, transducers: torch.Tensor | vedo.mesh.core.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, transducer_norms=None): View Source

272def BEM_forward_model_second_derivative_mixed(points:Tensor, scatterer:Mesh, transducers:Tensor|Mesh=None, use_cache_H:bool=True, k=Constants.k,
273                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False, 
274                           path:str="Media", p_ref=Constants.P_ref,internal_points = None, transducer_radius=Constants.radius, transducer_norms=None):
275    
276       
277    if transducers is None:
278        transducers = TRANSDUCERS
279
280    if H is None:
281        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, norms=transducer_norms)
282
283    Fxy, Fxz, Fyz = forward_model_second_derivative_mixed(points,transducers=transducers, k=k, transducer_radius=transducer_radius, p_ref=p_ref, transducer_norms=transducer_norms)
284    Gxy, Gxz, Gyz = get_G_second_mixed(points, scatterer, transducers, k=k)
285
286
287    Exy = Fxy + Gxy@H
288    Exz = Fxz + Gxz@H
289    Eyz = Fyz + Gyz@H
290
291    return Exy, Exz, Eyz

def BEM_laplacian( points: torch.Tensor, scatterer: vedo.mesh.core.Mesh, transducers: torch.Tensor | vedo.mesh.core.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, transducer_norms=None): View Source

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

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