In [1]:

import bempp.api
import numpy as np
import fractalpp

Set problem parameters:

In [2]:

prefractal_level = 3
#inc wave params:
kwave = 10
#mesh parameter
h = 2*np.pi/(10*kwave) # ten wavelengths long
inc_dir = np.sqrt(1/3)*np.array([1,1,-1])

Now create the mesh using fractalpp.mesh.koch We are including the optional parameter h_max, which will keep subdividing until all elements are less than this diameter. Other optional paramteres, and their defaults, are: h_max=math.inf, scale = 1, apex_angle = math.pi/3, shift = np.array([0,0]), mid_refine=True The final one of these avoivds refinement in the centre of the fractal, and only refines close to the edges.

In [3]:

koch_mesh = fractalpp.mesh.koch(prefractal_level,h_max = h)

Now plot the mesh. When prefractal_level is high, I would advise not running the below cell of code, as the plot function will takes ages and it's not necessary.

In [4]:

koch_mesh.plot()

Now make a function for plotting slices in this notebook, which we can use for each BC. We can also obtain our approximation to $u^s(x)\approx\mathcal{D}\phi(x) - \mathcal{S}\psi(x)$. The next few cells of code plot the solution in a slice $[-1,1]\times[-1,1]\times\{z\}\subset\mathbb{R}^3$ for some $z$.

In [5]:

def slice(u_s,z_val = -0.3):
    # z_val is the z value of square slice
    n_grid_points = 150
    plot_grid = np.mgrid[-1:1:n_grid_points*1j, -1:1:n_grid_points*1j]
    points = np.vstack((plot_grid[0].ravel(),
                        plot_grid[1].ravel(),
                        z_val*np.ones(plot_grid[0].size)))
    scattered_field_at_points = (u_s(points)).reshape((n_grid_points,n_grid_points))
    %matplotlib inline
    from matplotlib import pylab as plt

    plt.imshow(np.abs(scattered_field_at_points), extent=(-1,1,-1,1))
    plt.title('Slice of scattered field')
    plt.colorbar()

Sound-soft problem¶

In [ ]:

%%capture
# create incident wave
uinc = fractalpp.solve.u_i(kwave,inc_dir)
# solve impedance problem via BIE formulation
sol = fractalpp.solve.dirichlet(koch_mesh,uinc)

In [ ]:

sol.neumann_jump.plot()

In [ ]:

slice(sol.u_s)

Impedance problem¶

Now solve the impedance problem on the mesh, using our new impedance BIE.

In [ ]:

#impedance parameters:
lambda_plus = kwave*(1.5+1.5j)
lambda_minus = kwave*(1+1j)

In [ ]:

%%capture
# create incident wave
uinc = fractalpp.solve.u_i(kwave,inc_dir)
# solve impedance problem via BIE formulation
imp_sol = fractalpp.solve.impedance(koch_mesh,uinc,lambda_plus,lambda_minus)

The above fractalpp.solve.impedance has the optional parameter Hassen=False, which if set to True will implement the BIE from the Hassen et al paper, which is not well-posed for $\lambda_-+\lambda_+$=0.

From the imp_sol object we have created, we can access the components of the solution $\left(\phi,\psi\right) \approx \left([u],[\partial_nu]\right)$, via imp_sol.Dirichlet_jump and imp_sol.Neumann_jump. For example, we can run the following to get a surface plot:

In [ ]:

imp_sol.Dirichlet_jump.plot()

In [ ]:

imp_sol.Neumann_jump.plot()

In [ ]:

slice(imp_sol.u_s,z_val=0.1)