|
1 | 1 | import numpy as np |
| 2 | +import numba |
2 | 3 |
|
3 | | - |
4 | | -def laplace_smooth_data2(mesh,data,kappa=0.05,dt=1.0,iter_total=150): |
| 4 | +def laplace_smooth_data2_slow(mesh,data,kappa=0.05,dt=1.0,iter_total=150): |
5 | 5 | """ |
6 | 6 | Apply Laplacian smoothing to a scalar field on a mesh using explicit diffusion. |
| 7 | + |
| 8 | + DEPRECATED. Use the version without _slow in the name |
7 | 9 |
|
8 | 10 | This function performs iterative smoothing by updating each node's value |
9 | 11 | toward the average of its neighbors. It corresponds to forward Euler |
@@ -48,6 +50,38 @@ def laplace_smooth_data2(mesh,data,kappa=0.05,dt=1.0,iter_total=150): |
48 | 50 | return smooth_data |
49 | 51 |
|
50 | 52 |
|
| 53 | + |
| 54 | +def laplace_smooth_data2(mesh,data,kappa=0.05,dt=1.0,iter_total=150): |
| 55 | + """ |
| 56 | + Apply Laplacian smoothing to a scalar field on a mesh using explicit diffusion. |
| 57 | +
|
| 58 | + This function performs iterative smoothing by updating each node's value |
| 59 | + toward the average of its neighbors. It corresponds to forward Euler |
| 60 | + integration of the diffusion equation. |
| 61 | +
|
| 62 | + Parameters |
| 63 | + ---------- |
| 64 | + mesh : object |
| 65 | + A mesh object with a `nodes` array and a method `get_neighbor_nodes(ndx)` |
| 66 | + returning neighbor node indices for a given node. |
| 67 | + data : ndarray of shape (N,) |
| 68 | + Initial scalar field values at each mesh node. |
| 69 | + kappa : float, optional |
| 70 | + Diffusivity coefficient. Controls the rate of smoothing. Default is 0.05. |
| 71 | + dt : float, optional |
| 72 | + Timestep used for each smoothing iteration. Default is 1.0. |
| 73 | + iter_total : int, optional |
| 74 | + Number of smoothing iterations to perform. Default is 150. |
| 75 | +
|
| 76 | + Returns |
| 77 | + ------- |
| 78 | + smooth_data : ndarray of shape (N,) |
| 79 | + The smoothed scalar field after `iter_total` iterations. |
| 80 | + """ |
| 81 | + n_nodes = mesh.nodes.shape[0] |
| 82 | + neighbors = [np.array(mesh.get_neighbor_nodes(i), dtype=np.int32) for i in range(n_nodes)] |
| 83 | + |
| 84 | + return smooth_kernel_numba(data.astype(np.float64), neighbors, kappa, dt, iter_total) |
51 | 85 |
|
52 | 86 |
|
53 | 87 | def laplace_smooth_with_vel2(mesh,data,vel,kappa=0.05,dt=1.,iter_total=1): |
@@ -152,3 +186,32 @@ def laplace_smooth_with_vel3(mesh,nlayer,data,vel,kappa=0.05,dt=1.,iter_total=1) |
152 | 186 | iter += 1 |
153 | 187 | smooth_data = zz |
154 | 188 | return smooth_data |
| 189 | + |
| 190 | + |
| 191 | + |
| 192 | + |
| 193 | +@numba.njit |
| 194 | +def smooth_kernel_numba(data, neighbors, kappa, dt, iter_total): |
| 195 | + n_nodes = len(data) |
| 196 | + result = data.copy() |
| 197 | + for _ in range(iter_total): |
| 198 | + tmp = result.copy() |
| 199 | + for i in range(n_nodes): |
| 200 | + nbs = neighbors[i] |
| 201 | + if len(nbs) == 0: |
| 202 | + continue |
| 203 | + v_ave = 0.0 |
| 204 | + for j in nbs: |
| 205 | + v_ave += result[j] |
| 206 | + v_ave /= len(nbs) |
| 207 | + tmp[i] = result[i] + dt * kappa * (v_ave - result[i]) |
| 208 | + result = tmp |
| 209 | + return result |
| 210 | + |
| 211 | + |
| 212 | + |
| 213 | + |
| 214 | + |
| 215 | + |
| 216 | + |
| 217 | + |
0 commit comments