<h1>Results</h1>
-Using the above code (without the PML) and a forcing term of a Hertzian dipole
-at the center, we have generated the following solution wave. The two complex
-plane solutions are followed by the two solutions in the real plane. Furthermore,
-a significant resonance is observed, causing distorted images and necessitating
-a PML boundary condition. <br />
-TODO: add the images <br />
-The solution of the same problem with a PML of strength 2 of radii 8 and 10
-is shown below (in the same order or complex and real solutions). Clearly,
-the PML significantly reduces the distortion in our solution.<br />
-TODO: add the images<br />
-Additionally, an interface is added at y = 0, and by adjusting the surface
-conductivity value and the position of our dipole, we arrive at a standing
-wave.<br />
-TODO: add the images
+The solution is written to a .vtk file with four components. These are the real and imaginary parts of the $E_x$ and $E_y$ solution waves. With the current setup, the output should read
+
+@code
+Number of active cells: 4096
+Number of degrees of freedom: 16640
+Program ended with exit code: 0
+@endcode
+
+We first run the code without the interface and with the dipole centered at $(0,0)$. To demonstrate the effect of the PML, we first set the strength to 0. Change the following parameters in the .prm file:
+
+@code
+# is there an interface?
+ set interface boolean = false
+
+# position of the dipole
+ set dipole position = 0, 0
+
+# strength of the PML
+ set strength = 0
+@endcode
+
+The imaginary component of the E_x wave is shown below. We observe that in absence of the PML, there is a lot of distortion and resonance (the real parts will not be generated without a PML). Now, we increase the strength of the PML to 4.
+
+@code
+# strength of the PML
+ set strength = 4
+@endcode
+
+As we see, a much more coherent image is generated with an appropriate PML.
+
+<table width="80%" align="center">
+ <tr>
+ <td align="center">
+ <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML0_nointerface.png" alt="Visualization of the solution of step-81 with no interface and PML strength 0" height="500">
+ </td>
+ <td></td>
+<td align="center">
+ <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML4_nointerface.png" alt="Visualization of the solution of step-81 with no interface and PML strength 4" height="500">
+ </td>
+ </tr>
+</table>
+
+Now, let's generate a standing wave by adding an interface at the center. In order to observe this effect, we offset the center of the dipole to $(0, 0.8)$ and set the interface boolean to $true$
+
+@code
+# is there an interface?
+ set interface boolean = true
+
+# position of the dipole
+ set dipole position = 0, 0.8
+@endcode
+
+If the PML strength is set to 0, we see the following output (in order from left to right, the images represent the imaginary components of $E_x$ and $E_y$ and the real components of $E_x$ and $E_y$):
+<table width="80%" align="center">
+ <tr>
+ <td align="center">
+ <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML0_imagEx.png" alt="Visualization of the solution of step-81 with an interface and PML strength 0" height="300">
+ </td>
+ <td></td>
+<td align="center">
+ <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML0_imagEy.png" alt="Visualization of the solution of step-81 with an interface and PML strength 0" height="300">
+ </td>
+ <td></td>
+ <td align="center">
+ <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML0_realEx.png" alt="Visualization of the solution of step-81 with an interface and PML strength 0" height="300">
+ </td>
+ <td></td>
+<td align="center">
+ <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML0_realEy.png" alt="Visualization of the solution of step-81 with an interface and PML strength 0" height="300">
+ </td>
+ </tr>
+</table>
+
+Once again, we see significant distortion and resonance but the standing wave is visible. Now, we increase the PML strength to 4 and generate the following images (in order from left to right, the images represent the imaginary components of $E_x$ and $E_y$ and the real components of $E_x$ and $E_y$):
+<table width="80%" align="center">
+ <tr>
+ <td align="center">
+ <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML4_imagEx.png" alt="Visualization of the solution of step-81 with an interface and PML strength 4" height="300">
+ </td>
+ <td></td>
+<td align="center">
+ <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML4_imagEy.png" alt="Visualization of the solution of step-81 with an interface and PML strength 4" height="300">
+ </td>
+ <td></td>
+ <td align="center">
+ <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML4_realEx.png" alt="Visualization of the solution of step-81 with an interface and PML strength 4" height="300">
+ </td>
+ <td></td>
+<td align="center">
+ <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML4_realEy.png" alt="Visualization of the solution of step-81 with an interface and PML strength 4" height="300">
+ </td>
+ </tr>
+</table>
+
+As expected, the PML removes the distortion and resonance. The standing wave is also dissipating and getting absorbed within the PML, and as we increase the PML strength, the standing wave will dissipate more within the PML ring.
+
+
+
+
+
+
dipole_position = Point<dim>(0., 0.8);
add_parameter("dipole position",
dipole_position,
- "posititon of the dipole");
+ "position of the dipole");
dipole_orientation = Tensor<1, dim, double>{{0., 1.}};
add_parameter("dipole orientation",
unsigned int refinements;
unsigned int fe_order;
unsigned int quadrature_order;
- unsigned int n_outputs;
+ bool interface_bool;
void parse_parameters_callback();
void make_grid();
// @sect4{The Constructor}
// The Constructor simply consists specifications for the mesh
// and the order of the fnite elements. These are editable through
- // the .prm file.
+ // the .prm file. The interface_bool can be modified to remove the
+ // interface and the standing wave.
template <int dim>
Maxwell<dim>::Maxwell()
scaling = 20;
add_parameter("scaling", scaling, "scale of the hypercube geometry");
- refinements = 10;
+ refinements = 8;
add_parameter("refinements",
refinements,
"number of refinements of the geometry");
add_parameter("quadrature order",
quadrature_order,
"order of the quadrature");
+
+ interface_bool = true;
+ add_parameter("interface boolean",
+ interface_bool,
+ "is there an interface?");
}
fe = std::make_unique<FESystem<dim>>(FE_NedelecSZ<dim>(fe_order), 2);
}
- // Make the mesh for the domain and generate the triangulation required
+ // Make the mesh for the domain and generate the triangulation required.
+ // Additionally, there is an interface added here to visualize
+ // a standing wave. To generate a solution without any interface,
+ // comment out lines 455-459.
+
template <int dim>
void Maxwell<dim>::make_grid()
{
GridGenerator::hyper_cube(triangulation, -scaling, scaling);
triangulation.refine_global(refinements);
-
+
+ if (interface_bool){
for (auto &cell : triangulation.active_cell_iterators())
if (cell->center()[1] > 0.)
cell->set_material_id(1);
else
cell->set_material_id(2);
+ }
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
+
constraints.clear();
+
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+ VectorTools::project_boundary_values_curl_conforming_l2(
+ dof_handler,
+ 0, /* real part */
+ dealii::ZeroFunction<dim>(2 * dim),
+ 0, /* boundary id */
+ constraints);
+ VectorTools::project_boundary_values_curl_conforming_l2(
+ dof_handler,
+ dim, /* imaginary part */
+ dealii::ZeroFunction<dim>(2 * dim),
+ 0, /* boundary id */
+ constraints);
+
constraints.close();
+
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
const auto inner_radius =
perfectly_matched_layer.inner_radius;
+ if (id == 0){
+ auto mu_inv = parameters.mu_inv(position, id);
+ auto epsilon = parameters.epsilon(position, id);
- 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_tensor(position);
-
- mu_inv = mu_inv / d;
- epsilon = invert(A) * epsilon * invert(B);
- };
+ 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_tensor(position);
- const auto normal = fe_face_values.normal_vector(q_point);
+ mu_inv = mu_inv / d;
+ epsilon = invert(A) * epsilon * invert(B);
+ };
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- {
- 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);
+ const auto normal = fe_face_values.normal_vector(q_point);
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
- const auto phi_j =
- real_part.value(j, q_point) +
- 1.0i * imag_part.value(j, q_point);
- const auto phi_j_T =
- tangential_part(phi_j, normal) *
- fe_face_values.JxW(q_point);
-
- const auto prod = mu_inv * epsilon;
- const auto sqrt_prod = prod;
-
- const auto temp =
- -1.0i *
- scalar_product((sqrt_prod * phi_j_T), phi_i_T);
- cell_matrix(i, j) += temp.real();
- } /* j */
- } /* i */
- } /* q_point */
+ 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)
+ {
+ const auto phi_j =
+ real_part.value(j, q_point) +
+ 1.0i * imag_part.value(j, q_point);
+ const auto phi_j_T =
+ tangential_part(phi_j, normal) *
+ fe_face_values.JxW(q_point);
+
+ const auto prod = mu_inv * epsilon;
+ const auto sqrt_prod = prod;
+
+ const auto temp =
+ -1.0i *
+ scalar_product((sqrt_prod * phi_j_T), phi_i_T);
+ cell_matrix(i, j) += temp.real();
+ } /* j */
+ } /* i */
+ }
+ } /* q_point */
}
else
{
}
};
-
+ // We use a direct solver from the SparseDirectUMFPACK to solve the system
template <int dim>
void Maxwell<dim>::solve()
{
A_direct.vmult(solution, system_rhs);
}
-
+ // The output is writted into a vtk file with 4 components
template <int dim>
void Maxwell<dim>::output_results()
{
projects / dealii.git / commitdiff