and distortion in our solutions, we are implementing a Perfectly Matched Layer
(PML) in the scattering configuration.
-The concept of a PML was pioneered by Bérenger and it is is an indispensable tool
-for truncating unbounded domains for wave equations and often used in the
-numerical approximation of scattering problems. It is essentially a thin layer with
-modified material parameters placed near the boundary such that all outgoing
-electromagnetic waves decay exponentially with no “artificial” reflection due to
-truncation of the domain.
-
-Our PML is essentially a concentric circle with modified material coefficients
-($\varepsilon_r, \mu_r, \sigma$). It is located in a small region near the boundary
-$\partial\Omega$ and the transformation of the material coordinates is chosen to
-be a function of the radial distance $\rho$ from the origin $e_r$. The normal field
-$\nu$ of $\Sigma$ is orthogonal to the radial direction $e_r$, which makes
-$\mathbf{J}_a \equiv 0$ and $\mathbf{M}_a \equiv 0$ within the PML.
+The concept of a PML was pioneered by Bérenger @cite Berenger1994
+and it is is an indispensable tool for truncating unbounded domains for
+wave equations and often used in the numerical approximation of scattering
+problems. It is essentially a thin layer with modified material parameters
+placed near the boundary such that all outgoing electromagnetic waves decay
+exponentially with no “artificial” reflection due to truncation of the
+domain.
+
+Our PML is a concentric circle with modified material coefficients
+($\varepsilon_r, \mu_r, \sigma$). It is located in a small region near the
+boundary $\partial\Omega$ and the transformation of the material
+coordinates is chosen to be a function of the radial distance $\rho$ from
+the origin $e_r$. The normal field $\nu$ of $\Sigma$ is orthogonal to the
+radial direction $e_r$, which makes $\mathbf{J}_a \equiv 0$ and
+$\mathbf{M}_a \equiv 0$ within the PML.
@htmlonly
<p align="center">
<p> Solution with no interface, Dirichlet boundary conditions and PML strength 0.</p>
</td>
<td></td>
- <td align="center">
+ <td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-nointerface_noabs_PML0.png" alt="Visualization of the solution of step-81 with no interface, no absorbing boundary conditions and PML strength 0" height="210">
- <p> Solution with no interface, absorbing boundary conditions and PML strength 0.</p>
- </td>
- <td></td>
- <td align="center">
+
<p> Solution with no interface, absorbing boundary conditions and PML strength 0.</p>
+ </td>
+ <td></td>
+ <td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-nointerface_abs_PML4.png" alt="Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4" height="210">
<p> Solution with no interface, absorbing boundary conditions and PML strength 4.</p>
</td>
<table width="80%" align="center">
<tr>
- <td align="center">
+ <td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-imagEx_noabs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
<p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
- </td>
- <td></td>
- <td align="center">
+ </td>
+ <td></td>
+ <td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-imagEx_abs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
- <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
+
<p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
</td>
<td></td>
<td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-imagEx_abs_PML4.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
- <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
+
<p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
</td>
</tr>
</table>
<table width="80%" align="center">
<tr>
- <td align="center">
+ <td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-realEx_noabs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
- <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
- </td>
- <td></td>
- <td align="center">
+
<p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
+ </td>
+ <td></td>
+ <td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-realEx_abs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
- <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
+
<p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
</td>
<td></td>
<td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-realEx_abs_PML4.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
- <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
+
<p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
</td>
</tr>
</table>
Here are some animations to demonstrate the effect of the PML
<table width="80%" align="center">
<tr>
- <td align="center">
+ <td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-dirichlet_Ex.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
<p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
- </td>
- <td></td>
- <td align="center">
+ </td>
+ <td></td>
+ <td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-absorbing_Ex.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
- <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
+
<p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
</td>
<td></td>
<td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-perfectly_matched_layer_Ex.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
- <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
+
<p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
</td>
</tr>
</table>
<table width="80%" align="center">
<tr>
- <td align="center">
+ <td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-dirichlet_Ey.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
<p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
- </td>
- <td></td>
- <td align="center">
+ </td>
+ <td></td>
+ <td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-absorbing_Ey.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
- <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
+
<p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
</td>
<td></td>
<td align="center">
<img src="https://www.dealii.org/images/steps/developer/step-81-perfectly_matched_layer_Ey.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
- <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
+
<p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
</td>
</tr>
</table>
// The set of include files is quite standard. The most notable include is
// the fe/fe_nedelec_sz.h file which allows us to use the FE_NedelecSZ elements.
// This is an implementation of the $H^{curl}$ conforming Nédélec Elements
-// that resolves the sign conflict issues that arise from parametrization.
+// that resolves the sign conflict issues that arise from parametrization
+// (for details we refer to the documentation of the FE_NedelecSZ element).
#include <deal.II/base/function.h>
#include <deal.II/base/parameter_acceptor.h>
class PerfectlyMatchedLayer : public ParameterAcceptor
{
public:
- static_assert(dim == 2, "dim == 2"); /* only works in 2D */
+ static_assert(dim == 2,
+ "The perfectly matched layer is only implemented in 2D.");
Parameters<dim> parameters;
PerfectlyMatchedLayer();
- double inner_radius;
- double outer_radius;
- double strength;
-
std::complex<double> d(const Point<dim> point);
std::complex<double> d_bar(const Point<dim> point);
rank2_type b_matrix(const Point<dim> point);
rank2_type c_matrix(const Point<dim> point);
+
+ private:
+ double inner_radius;
+ double outer_radius;
+ double strength;
};
typename std::complex<double>
PerfectlyMatchedLayer<dim>::d(const Point<dim> point)
{
- const auto radius = point.norm();
- const double s =
- strength * ((radius - inner_radius) * (radius - inner_radius)) /
- ((outer_radius - inner_radius) * (outer_radius - inner_radius));
- return 1.0 + 1.0i * s;
+ const auto radius = point.norm();
+ if (radius > inner_radius)
+ {
+ const double s =
+ strength * ((radius - inner_radius) * (radius - inner_radius)) /
+ ((outer_radius - inner_radius) * (outer_radius - inner_radius));
+ return 1.0 + 1.0i * s;
+ }
+ else
+ {
+ return 1.0;
+ }
}
typename std::complex<double>
PerfectlyMatchedLayer<dim>::d_bar(const Point<dim> point)
{
- const auto radius = point.norm();
- const double s_bar =
- strength / 3. *
- ((radius - inner_radius) * (radius - inner_radius) *
- (radius - inner_radius)) /
- (radius * (outer_radius - inner_radius) * (outer_radius - inner_radius));
- return 1.0 + 1.0i * s_bar;
+ const auto radius = point.norm();
+ if (radius > inner_radius)
+ {
+ const double s_bar =
+ strength / 3. *
+ ((radius - inner_radius) * (radius - inner_radius) *
+ (radius - inner_radius)) /
+ (radius * (outer_radius - inner_radius) *
+ (outer_radius - inner_radius));
+ return 1.0 + 1.0i * s_bar;
+ }
+ else
+ {
+ return 1.0;
+ }
}
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const Point<dim> &position = quadrature_points[q_point];
- const auto radius = position.norm();
- const auto inner_radius = perfectly_matched_layer.inner_radius;
auto mu_inv = parameters.mu_inv(position, id);
auto epsilon = parameters.epsilon(position, id);
const auto J_a = parameters.J_a(position, id);
- if (radius >= inner_radius)
- {
- auto A = perfectly_matched_layer.a_matrix(position);
- auto B = perfectly_matched_layer.b_matrix(position);
- auto d = perfectly_matched_layer.d(position);
+ const auto A = perfectly_matched_layer.a_matrix(position);
+ const auto B = perfectly_matched_layer.b_matrix(position);
+ const auto d = perfectly_matched_layer.d(position);
- mu_inv = mu_inv / d;
- epsilon = invert(A) * epsilon * invert(B);
- };
+ mu_inv = mu_inv / d;
+ epsilon = invert(A) * epsilon * invert(B);
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (const auto i : fe_values.dof_indices())
{
const auto phi_i = real_part.value(i, q_point) -
1.0i * imag_part.value(i, q_point);
(1.0i * scalar_product(J_a, phi_i)) * fe_values.JxW(q_point);
cell_rhs(i) += rhs_value.real();
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ for (const auto j : fe_values.dof_indices())
{
const auto phi_j = real_part.value(j, q_point) +
1.0i * imag_part.value(j, q_point);
++q_point)
{
const auto &position = quadrature_points[q_point];
- const auto radius = position.norm();
- const auto inner_radius =
- perfectly_matched_layer.inner_radius;
auto mu_inv = parameters.mu_inv(position, id);
auto epsilon = parameters.epsilon(position, id);
- if (radius >= inner_radius)
- {
- auto A = perfectly_matched_layer.a_matrix(position);
- auto B = perfectly_matched_layer.b_matrix(position);
- auto d = perfectly_matched_layer.d(position);
+ const auto A =
+ perfectly_matched_layer.a_matrix(position);
+ const auto B =
+ perfectly_matched_layer.b_matrix(position);
+ const auto d = perfectly_matched_layer.d(position);
- mu_inv = mu_inv / d;
- epsilon = invert(A) * epsilon * invert(B);
- };
+ mu_inv = mu_inv / d;
+ epsilon = invert(A) * epsilon * invert(B);
const auto normal =
fe_face_values.normal_vector(q_point);
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (const auto i : fe_face_values.dof_indices())
{
const auto phi_i =
real_part.value(i, q_point) -
1.0i * imag_part.value(i, q_point);
const auto phi_i_T = tangential_part(phi_i, normal);
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ for (const auto j : fe_face_values.dof_indices())
{
const auto phi_j =
real_part.value(j, q_point) +
++q_point)
{
const auto &position = quadrature_points[q_point];
- const auto radius = position.norm();
- const auto inner_radius =
- perfectly_matched_layer.inner_radius;
auto sigma = parameters.sigma(position, id1, id2);
- if (radius >= inner_radius)
- {
- auto B = perfectly_matched_layer.b_matrix(position);
- auto C = perfectly_matched_layer.c_matrix(position);
- sigma = invert(C) * sigma * invert(B);
- };
+ const auto B = perfectly_matched_layer.b_matrix(position);
+ const auto C = perfectly_matched_layer.c_matrix(position);
+ sigma = invert(C) * sigma * invert(B);
const auto normal = fe_face_values.normal_vector(q_point);
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (const auto i : fe_face_values.dof_indices())
{
const auto phi_i = real_part.value(i, q_point) -
1.0i * imag_part.value(i, q_point);
const auto phi_i_T = tangential_part(phi_i, normal);
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ for (const auto j : fe_face_values.dof_indices())
{
const auto phi_j =
real_part.value(j, q_point) +
projects / dealii.git / commitdiff