<jupyter_start><jupyter_text>Software profesional en Acústica 2020-21 (M2i)*This notebook contains a modification of the notebook [FEM_Helmholtz_equation_Robin](https://github.com/spatialaudio/computational_acoustics/blob/master/FEM_Helmholtz_equation_Robin.ipynb), created by <NAME>, <NAME>, Computational Acoustics Examples, 2018. The text/images are licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/). The code is released under the [MIT license](https://opensource.org/licenses/MIT).*  Numerical Solution of the Helmholtz Equation with Robin Boundary Conditions using the Finite Element Method This notebook illustrates the numerical solution of the wave equation for an harmonic excitation and Robin boundary conditions using the [Finite Element Method](https://en.wikipedia.org/wiki/Finite_element_method) (FEM). The method aims at an approximate solution by subdividing the area of interest into smaller parts with simpler geometry, linking these parts together and applying methods from the calculus of variations to solve the problem numerically. The FEM is a well established method for the numerical approximation of the solution of partial differential equations (PDEs). The solutions of PDEs are often known analytically only for rather simple geometries. FEM based simulations allow to gain insights into other more complex cases.  Problem StatementThe inhomogeneous [Helmholtz equation](https://en.wikipedia.org/wiki/Helmholtz_equation) is given as\begin{equation}\Delta P(\mathbf{x}, \omega) + \frac{\omega^2}{c^2} P(\mathbf{x}, \omega) = - Q(\mathbf{x}, \omega) .\end{equation}We aim for a numerical solution of the Helmholtz equation on the domain $V$ with respect to the homogeneous Robin boundary condition\begin{equation}V_n(\mathbf{x}, \omega) + \frac{1}{Z(\mathbf{x}, \omega)} P(\mathbf{x}, \omega) = 0 \qquad \text{for } x \in \partial V ,\end{equation}where $V_n(\mathbf{x}, \omega)$ denotes the particle velocity in inward normal direction to the boundary $\partial V$ of $V$ and $Z(\mathbf{x}, \omega)$ the acoustic impedance of the boundary.The particle velocity can be linked to the pressure using the Euler equation\begin{equation}-\mathrm{j} \omega \rho_0 V_n(\mathbf{x}, \omega) = \frac{\partial}{\partial n} P(\mathbf{x}, \omega) ,\end{equation}where $\rho_0$ denotes the static density of air.Introducing this into the Robin boundary equation above yields\begin{equation}\frac{\partial}{\partial n} P(\mathbf{x}, \omega) - \mathrm{j} \underbrace{\frac{\omega \rho_0}{Z}}_{\sigma} P(\mathbf{x}, \omega) = 0 \qquad \text{for } x \in \partial V .\end{equation}The medium impedance of air for free-field propagation is $Z_0 = \rho_0 c$, hence $\sigma_0 = \frac{\omega}{c}$ in this case. Free-field conditions can be simulated by matching the impedance of the boundary to $Z_0$.  Variational FormulationStarting from the [variational formulation of the Helmholtz equation](FEM_Helmholtz_equation_2D.ipynbVariational-Formulation) (before application of Green's first theorem)\begin{equation}{-} \int_V \nabla P(\mathbf{x}, \omega) \cdot \nabla V(\mathbf{x}, \omega) \mathrm{d}x  + \int_{\partial V} V(\mathbf{x}, \omega) \frac{\partial}{\partial n}  P(\mathbf{x}, \omega) \mathrm{d}s + \frac{\omega^2}{c^2} \int_V P(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x = -\int_V Q(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x\end{equation}and introducing the Robin boundary condition into the second integral yields\begin{equation}{-} \int_V \nabla P(\mathbf{x}, \omega) \cdot \nabla V(\mathbf{x}, \omega) \mathrm{d}x  + \mathrm{j} \sigma \int_{\partial V} V(\mathbf{x}, \omega) P(\mathbf{x}, \omega) \mathrm{d}s + \frac{\omega^2}{c^2} \int_V P(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x = -\int_V Q(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x .\end{equation}It is common to express this integral equation in terms of the bilinear $a(P, V)$ and linear $L(V)$ forms \begin{equation}a(P, V) = \frac{\omega^2}{c^2} \int_V P(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x - \int_V \nabla P(\mathbf{x}, \omega) \cdot \nabla V(\mathbf{x}, \omega) \mathrm{d}x + \mathrm{j} \sigma \int_{\partial V} V(\mathbf{x}, \omega) P(\mathbf{x}, \omega) \mathrm{d}s , \end{equation}\begin{equation}L(V) = -\int_V Q(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x ,\end{equation}where\begin{equation}a(P, V) = L(V) .\end{equation}Computational implementations of the FEM (like FEniCS) may not be able to handle complex numbers.In this case the problem can be split into two coupled problems with respect to the real and imaginary part. By introducing $P(\mathbf{x}, \omega) = P_r(\mathbf{x}, \omega) + \mathrm{j} P_i(\mathbf{x}, \omega)$ and $V(\mathbf{x}, \omega) = V_r(\mathbf{x}, \omega) + \mathrm{j} V_i(\mathbf{x}, \omega)$ and identifying the real and imaginary parts of the bilinear and linear forms, we get\begin{equation}a_r = \int_V \left( \frac{\omega^2}{c^2} V_r(\mathbf{x}, \omega) P_r(\mathbf{x}, \omega) -  \frac{\omega^2}{c^2} V_i(\mathbf{x}, \omega) P_i(\mathbf{x}, \omega) - \nabla P_r(\mathbf{x}, \omega) \cdot \nabla V_r(\mathbf{x}, \omega) + \nabla P_i(\mathbf{x}, \omega) \cdot \nabla V_i(\mathbf{x}, \omega) \right) \mathrm{d}x -\sigma \cdot \int_{\partial V} \left( V_i(\mathbf{x}, \omega) P_r(\mathbf{x}, \omega) + V_r(\mathbf{x}, \omega) P_i(\mathbf{x}, \omega) \right) \mathrm{d}s, \end{equation}\begin{equation}a_i = \int_V \left( \frac{\omega^2}{c^2} V_r(\mathbf{x}, \omega) P_i(\mathbf{x}, \omega) +  \frac{\omega^2}{c^2} V_i(\mathbf{x}, \omega) P_r(\mathbf{x}, \omega) - \nabla P_i(\mathbf{x}, \omega) \cdot \nabla V_r(\mathbf{x}, \omega) + \nabla P_r(\mathbf{x}, \omega) \cdot \nabla V_i(\mathbf{x}, \omega) \right) \mathrm{d}x +\sigma \cdot \int_{\partial V} \left( V_r(\mathbf{x}, \omega) P_r(\mathbf{x}, \omega) - V_i(\mathbf{x}, \omega) P_i(\mathbf{x}, \omega) \right) \mathrm{d}s\end{equation}for the bilinear form.  Numerical SolutionThe numerical solution of the variational problem is based on [FEniCS](https://fenicsproject.org/), an open-source framework for numerical solution of PDEs.Its high-level Python interface `dolfin` is used in the following to define the problem and compute the solution.The implementation is based on the variational formulation derived above.It is common in the FEM to denote the solution of the problem by $u$ and the test function by $v$.The definition of the problem in FEniCS is very close to the mathematical formulation of the problem.For the subsequent examples the solution of inhomogeneous wave equation for a point source $Q(\mathbf{x}) = \delta(\mathbf{x}-\mathbf{x_s})$ at position $\mathbf{x_s}$ is computed using the FEM.A function is defined for this purpose, as well as for the plotting of the resulting sound field.<jupyter_code>import dolfin
import mshr
import matplotlib.pyplot as plt
%matplotlib inline


def Helmholtz_Robin(mesh, frequency, xs, sigma=dolfin.Constant(0), c=343):
    
    # squared wavenumber
    k2 = dolfin.Constant(2*dolfin.pi*frequency/c)**2
    
    # define function space
    V = dolfin.VectorFunctionSpace(mesh, "CG", 1, dim=2)

    # define variational problem
    (u_r, u_i) = dolfin.TrialFunction(V)
    (v_r, v_i) = dolfin.TestFunction(V)

    a_r = ( k2 * dolfin.inner(u_r,v_r) - k2 * dolfin.inner(u_i,v_i) - dolfin.inner(dolfin.grad(u_r), dolfin.grad(v_r)) + dolfin.inner(dolfin.grad(u_i), dolfin.grad(v_i)) ) * dolfin.dx  - sigma*dolfin.inner(u_r, v_i) * dolfin.ds - sigma*dolfin.inner(u_i, v_r) * dolfin.ds
    a_i = ( k2 * dolfin.inner(u_r,v_i) + k2 * dolfin.inner(u_i,v_r) - dolfin.inner(dolfin.grad(u_r), dolfin.grad(v_i)) - dolfin.inner(dolfin.grad(u_i), dolfin.grad(v_r)) ) * dolfin.dx  + sigma*dolfin.inner(u_r, v_r) * dolfin.ds - sigma*dolfin.inner(u_i, v_i) * dolfin.ds
    L_r = dolfin.Constant(0) * v_r * dolfin.dx
    L_i = dolfin.Constant(0) * v_i * dolfin.dx

    a = a_r + a_i
    L = L_r + L_i
    A, b = dolfin.assemble_system(a, L)

    # define inhomogenity
    delta = dolfin.PointSource(V.sub(0), xs, -1)  # negative amplitude accounts for -Q(x,w) in inhomogeneous wave equation
    delta.apply(b)

    # compute solution
    u = dolfin.Function(V)
    dolfin.solve(A, u.vector(), b)
    (u_r, u_i) = dolfin.split(u)
    
    return u_r


def plot_soundfield(u):
    '''plots solution of FEM-based simulation'''
    fig = plt.figure(figsize=(10,10))
    fig = dolfin.plot(u)
    plt.title(r'$P(\mathbf{x}, \omega)$')
    plt.xlabel(r'$x$ in m')
    plt.ylabel(r'$y$ in m')
    plt.colorbar(fig, fraction=0.038, pad=0.04);<jupyter_output><empty_output><jupyter_text>Sound Field in a Rectangular RoomThe two-dimensional sound field in a rectangular room (rectangular plate) with homogeneous Robin boundary conditions is computed for the frequency $f=120$ Hz and source position $x_s = (1.2,3.2)$ m.<jupyter_code>f = 120  # frequency
xs = dolfin.Point(1.2, 3.2)  # source position

# define geometry and mesh
mesh = dolfin.RectangleMesh(dolfin.Point(0, 0), dolfin.Point(5, 4), 200, 200, "right/left")<jupyter_output><empty_output><jupyter_text>First, the case of sound-hard (Neumann) boundary conditions with $\sigma = 0$ ($Z \to \infty$) is considered.<jupyter_code># compute solution for sigma=0
u = Helmholtz_Robin(mesh, f, xs, sigma=dolfin.Constant(0))

# plot sound field
plot_soundfield(u)
plot_soundfield(abs(u))
plt.title(r'$|P(\mathbf{x}, \omega)|$');<jupyter_output>Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.<jupyter_text>Now the case of a matched boundary (free-field propagation) with $\sigma = \frac{\omega}{c}$ ($Z = Z_0$) is considered.<jupyter_code># compute solution for free-field propogation
u = Helmholtz_Robin(mesh, f, xs, sigma=dolfin.Constant(2*dolfin.pi*f/343))

# plot sound field
plot_soundfield(u)
plot_soundfield(abs(u))
plt.title(r'$|P(\mathbf{x}, \omega)|$');<jupyter_output><empty_output><jupyter_text>The last example shows the simulation result for sound-soft (Dirichlet) boundary conditions with $\sigma \to \infty$ ($Z=0$).<jupyter_code># compute solution for very large sigma
u = Helmholtz_Robin(mesh, f, xs, sigma=dolfin.Constant(1e15))

# plot sound field
plot_soundfield(u)
plot_soundfield(abs(u))
plt.title(r'$|P(\mathbf{x}, \omega)|$');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1529569122002484, 'fraction_numerical': 0.010127066017005828, 'mean_word_length': 4.7548103353490925, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '624859', 'license': 'MIT', 'n_tokens_mistral': 3708, 'n_tokens_neox': 3601, 'n_words': 1062, 'path': 'notebooks/FEniCS_Helmholtz_equation_Robin_complex.ipynb', 'repo_name': 'maprieto/SoftwareProfesionalAcustica', 'size': '10480'}